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Summary 

The computerized form of the minimum 
commitment method of interpolating and 
extrapolating stress versus time to failure data is a 
program called MEGAI6. This report describes 
MEGA 16, giving examples of its many plots and 
tabular outputs for a typical set of data. The 
program assumes a specific model equation and then 
provides a family of predicted isothermals for any set 
of data with at least 12 stress-rupture results from 
three different temperatures spread over reasonable 
stress and time ranges. It is written in Fortran IV 
using IBM plotting subroutines, and it runs on an 
IBM 370 time sharing system. 


Introduction 

Fundamental to current design practice of ground 
power turbine equipment is an estimate of a 
maximum allowable stress at which a particular 
material can be expected to survive for at least 
100 000 hours (11.6 yr) of exposure to a specified 
temperature. The search for a simple, reliable 
method of predicting such an estimated stress has 
been pursued for many years, by many individuals, 
using many scientific and empirical techniques (refs. 
1 and 2). One of the more common approaches 
involves the use of a time-temperature parameter to 
extrapolate data from the standard stress-rupture test 
results. In such a method a model equation relating 
observed rupture life to temperature and stress is 
assumed, and its unknown terms are estimated using 
available test data. Characteristically, the data for a 
particular material are sparse, cover limited ranges of 
stress and temperature, or else are collected from a 
mixture of different heats, product forms, or other 
manufacturing and testing variables. In addition, the 
model itself introduces other assumptions and 
oversimplifications. Nonetheless, a time-temperature 
parameter is expected to furnish a reasonable 
estimate of maximum allowable stress at a given life 
or life at a given stress. This can be used in the design 
of machine components that are required to operate 
at high temperatures for extended periods. 

The purpose of this paper is to present the 
computer program, MEGAI6 (Manson Ensign 
Generalized Analysis, Version 16), which has been 
developed as an objective, fast, and reliable way to 
use a particular time-temperature method. The 


method, often referred to as the “minimum 
commitment method,” has been described previously 
(ref. 3). Results from its application to many sets of 
data using MEGAI6 have been published (ref. 4); they 
compare favorably with results from other manual 
and computerized methods. A complete listing of 
MEGA 16 and its subroutines is given in the appendix. 


Program Description 


General Features 

Basically, mega 1 6 is a computer program which 
fits, to one or many sets of input data, a model 
equation of the form 

log t+A(9 log /-h(P = 8 

where (P and 8 are temperature and stress functions, 
each with several coefficients to be determined. 
MEGAI 6 produces a series of plots that represent these 
functions and the stress rupture behavior of the 
material over the range of input data and to an 
extrapolated time to failure of up to 100 (XX) hours. 

There are three sections to the main program (fig. 1): 
The first contains all the general housekeeping 
(dimensions, constants, etc.) and reading of data and 
variables associated with a given set of data. The 
second has the main loop which sets up and solves the 
equations for each value of an adjustable parameter 
called A, calculates the predicted values and 
statistics, prints the tabular results, and plots the 
curves for that A value. The third section produces 
other ancillary plots to show the merits of various A 
values used and the residuals of logarithm of time to 
failure versus other variables to indicate any trends 
that may have resulted with the model equation and 
the given data. 

After looping through the various A values 
selected in the second section of the program (plus an 
extra A value read in the data) and doing the final 
plots in the third section, the program then goes on to 
the next set of data for another material and repeats 
the whole procedure. Thus a great number of data 
sets may be handled with one loading of the 
program. 

Published stress-rupture data usually include fewer 
than 100 specimens for a material. The present form 
of MEGAI 6 can handle 200 data points per set; this 
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cards and raw data. It also affords considerable 
room for expansion or addition of other 
information, such as creep behavior, tensile strength, 
or even comments about a particular data point. The 
first nine lines of every data set give the title, 
Section I counters, constants, plotting values, and the number 

5 Do loops follow. Each datum is then given on 

succeeding lines, with (from left to right) the values 
for temperature (in deg F), stress (in ksi), hours to 
failure, and number of times the point is to be 
weighted, along with other comments or information 
that will not be read by the program. The data are 
usually listed in the order of increasing temperature 
and decreasing stress. Also, all data with the same 
temperature are grouped together, while data with 
the same temperature and stress value are grouped in 
order of increasing life. Of course these formats 
could be changed to match a given set of data, but we 
Section II have found it preferable to use a second program 

31 Do loops which transforms data from another format to this 

format. At the present, many data sets, for over 30 
different materials, are stored on tape in this form, 
ready to run. 

Calculations 


Yes 



The general equation used for mega 1 6 (ref. 4) is 

Plot A scan and residuals 



log t+A 

log t (P(7) + (P(7) = g(log S) (1) 



y Section III 
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9(logS) 

a stress function 


Figure 1. - Flowchart for MEGA16. 


S a/Of^ 

o a given stress value 

Oq a reference stress, near midrange 


value could be increased, if necessary, by appropriate 
expansion of dimensions. The maximum size of the 
matrix to be inverted is 200 by 48. For this model 
equation the 48 dimension could be reduced 
considerably since there are not as many terms as in 
other model equations using this same framework. 
Since mega 1 6 has been developed on a computer 
which has virtual memory, minimal effort was 
devoted to minimize storage requirements. 

Input 


The constant A takes on selected values from -0.2 
to 0.2 and also can be modified by the following 
equation: 


A*=A 


1 - AKON 


^ '^mid 


2 


( 2 ) 


where 

T a given temperature value 

a reference temperature at midrange 
AKON a constant, usually 0, 15, or 30 


Figure 2 is an example, with annotations, of the 
standard input data format required by MEGA16. The 
format is not efficient in terms of storage space 
required, but is has been very convenient for sorting, 
editing, and changing the constants of the control 


The temperature function (P(7) is expressed in the 
form of a station function, which is mathematically 
similar to a LaGrangian interpolation (see fig. 3 and 
ref. 3). The stress function is composed of two parts, 
each of the form 
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(VARIABLE NAMES 


FORMAT 


33 ASTROLOY MEGA14 TITEL 16A4 


0 

3 

0 1 





NF NP NG 

MK 

1615 

5 

5 

3 5 

1 




NL NT NS 

NISO NXAV 

1615 

1400 

1500 

1600 1700 

1800 




Cl SOT 


16F5.0 

1400 

1500 

1600 1700 

1800 




ALOGT=TSTA 


16F5.0 

1400 

1500 

1600 1700 

1800 




TSTA 


16F5.0 

1.16 

1.632 

• 004 





SSTA 


16F5.0 

33 

33 

99000 - 

-0.15 


FTFTFTFTTFTTTFTTT 

NODU NODT 

CUTOF AEXT 

OPT 2I5»F9.1»F6.3»4X»17L1 

100 

80 

60 50 

40 


30 20 10 

1.0 -1.0 

CONSTR <rields 9*10 = ALPHA) 16F5.0 

1400101. 

12.8 

0 




TMP SIG 

AT UF BATCH 

F5 . O » F7 . 3 f F8 ♦ 2 * 2F5 . O 

1400 

B6. 

59. 

0 







1400 

80. 

176.6 

0 







1400 

74. 

400.7 

0 







1400 

70. 

577. 

0 







1400 

61. 

2279.8 

0 







1400 

55. 

4063.2 

0 







1500 

75. 

30.5 

0 







1500 

64. 

142.2 

0 







1500 

56. 

351.3 

0 




COLUMN 

VARIABLE 

UNITS 

1500 

52. 

712. 

0 







1500 

45. 

1228.3 

0 




1 

t,eMF>er3t.ure 

decrees F 

1500 

3^. 

2227 . 4 

0 




2 

s-tress 

KSI 


31. 

4393.4 

0 




3 

Liee 

hours 

1600 

64. 

10.5 

0 




4 

weidhLind 


1600 

56.5 

28.8 

0 







1600 

46.5 

145.8 

0 


► t.o^3l or NODU 

values 




1600 

41. 

253.0 

0 







1600 

37. 

535.7 

0 







1600 

31. 

888. 

0 







1600 

24.5 

2899.7 

0 







1600 

19. 

6331 . 

0 







1700 

41 . 

11.5 

0 







1700 

33.5 

44.2 

0 







1700 

29. 

120.9 

0 







1700 

24. 

342.7 

0 







1700 

21. 

746.7 

0 







1700 

17.5 

1768.7 

0 







1700 

14.5 

2838 . 7 

0 







laoo 

29.5 

6.1 

0 







leoo 

20.5 

49.3 

0 







leoo 

17. 

174. 

0 







1800 

14.5 

340.7 

0 








Figure 2 . - Typical data set and input format for MEGA16, 


9(log S) = £) + £log 5 + 
where 

a an integer, usually + 

- 1 for (7>(Tq 

The two parts are connected at a spline point where 
the second derivatives with respect to log S are forced 
to be equal and will also be continuous if -o;. 

The calculations of megai6 begin with the 
expression of each of the input data points in this 
form of the model equation 


(3) 


1 for o<Op and 


log // = (P(7)[/4* log^ + l] + £>+£'log(^)+Fi 


(g;/go)“ _ 1 _ log(g//qo) ' 

. 5.30«2 5.30a2 2.30a 


p f (q//go)“ _ 1 _ lo g(g//go) ' 

5.30a2 5.30a2 2.30(-a) . 


(4) 


where 

F\ is set equal to 0 for ^><7^ 
F 2 is set equal to 0 a<ao 


For the number of data used, / = NODU, and the 
model equation, a matrix of NODU by (NPAR + 1) 
dimensions is formed and then inverted to provide 
estimates for the coefficients D, E, F\, and F 2 , as 
well as the station function values representing the NP 
discrete temperatures. Figure 4 shows, for a set of 
data with unevenly spaced test temperatures, an 
example of this input matrix. The solution, by 
MEGA 16, of this set of equations representing the 
iron-nickel alloy A-286 (ref. 1) gave these 
coefficients: /?/=-2.17, -0.81, 0.0, 1.39, 2.64; 
D = 3.19; £’=-4.33; £1 =2.25; and £2 =-61.41. 
Predicted values of the logarithm of time to rupture 
are made by MEGAI6 using the coefficients with 
combinations of temperature and stress increments. 
To provide the predicted isothermals, equation (4) is 
rewritten as 


lest ~~ Qest ” 1 " ^est 


(5) 


where 


^est =£1 T’l +P 2 T 2 + ^PnpTnp (6) 


3 



P(Tt * Pn_i + Pn K„ + Pn+j + p„+2 l<n+2 


. I <r-V(T-Tn,i) 


K 


n 


1 

2 


K 


. 1 
n+1 ■ 2 


T-To-lW-T n.lWn ‘ ^0.2' * ^ ‘ 'n-l' 

•^n •'n-l"'n •fnn'ffn '’nja' 
n-Tn-il(T-Tn l(Tn,i - T„,2l • ff - T„I(I - I„»2«T„,i - T„.,l 
•"n • To*l'*'n*l ■ Tn-l^n*! ' V»2' 



Qest =/^+£’ + log 


(?) 


I p.\ (oi/<yo)°‘ _ 1 _ log(q//go) 

L 5.30a2 5.30a2 2.30a 


(7) 


Standard deviation = 


NODU 

(log t, -log ti)^ 

/=! 

NODU - NPAR - I 


(8) 


where 

log tj logarithm of observed life values 

log f| logarithm of predicted life values 

NODU number of data used 

NPAR number of parameters 

Output 

Figure 5 gives examples of the long and short 
tabular outputs from a run using the input data from 
figure 2. Figures 6 shows some of the 16 frames that 
were plotted on microfiche using an adaptation of an 
IBM system of plotting (ref. 5). Five plots for each of 
three A values plus the one A scan were requested by 
the series of F*s and T’s (F gives a plot, T omits it) in 
the eighth line of input for this material (fig. 2). The 
first two A values (0 and -0.5) were specified by a 
statement in the main program of megai6, and the 
third /!( - 0. 1 5) was specified in the data set (also line 
8 of fig. 2). 

Experience has shown that analysis of a given 
material is seldom, if ever, completed by the first run 
of that data set. This, coupled with the fact that 
constant A in the model equation of megai6 needs to 
be optimized, leads to the dual system of selecting 
trial A values; that is, a few values over the range of 
possible A values (-0.2 to 0.2) are set in the 
beginning of the main program as starting values; 0.1 
and -0.05 are useful for a new set of data. A plot of 
the standard deviation values versus a few A values 
often indicates a U-shaped pattern, pointing to a 
possible minimum standard deviation at another A 
value. This other A value can then be the one read in 
with the data as the “extra ^4” value (line 8 of fig. 2 
and -0.15 above) for succeeding runs. The short 
form of the tabular output lists only the standard 
deviation results for specified A values and thus is 
used without any other output as a quick means of 
finding the best A value. 


This permits log test values to be calculated for 
evenly spaced values of log(a//ao), which are then 
plotted versus log stress. 

Another important part of the calculation is the 
estimation of error. The standard deviation of the 
data used gives an idea of how well the model 
equation fits over the range of the input variables. It 
is calculated using 


Computer Program Details 

MEGA16 consists of a main program of 1300 lines 
and subroutines totaling 500 lines, including many 
for documentation. For a typical single set of data on 
an IBM 370 time sharing system, 15 seconds of 
central processing time are required, but this drops 
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Matrix column number 1 2345678 9 

Row T, 0 , t, logt *^5 I 9l I 92 I 93 Tg^ 

ksi hr 

1 1050 76 81.4 1.911 0.904 ~0 To 0 1.0 0.32lTo 0.041 

2 72 97. 9 1. 991 . 900 . 297 . 036 

3 69 210. 2 2. 323 . 834 . 279 . 032 

4 67 523. 3 2. 719 . 864 . 266 . 029 

5 56 12833.6 4.108 .795 | .188 .015 

6 1100 68 35^1 2.547 . 291 0.873 . 273 . 030 

7 1150 65 41.9 1.622 0 . 919 . 253 . 027 

8 62 177. 7 2. 250 . 888 . 232 . 023 

9 54 644.4 ^ 809 . 860 .172 . 013 

10 43 675^ 8 3.829 . 808 . 073 . 003 

11 38 15460.8 4.189 . 791 .020 . 0002 

12 1200 60 18.9 1.276 0 .218 .020 

13 56 154.1 Zm .188 .015 

14 54 385.6 i586 .172 .013 

15 47 839.0 2.924 .112 .006 

16 35 6882.4 3.838 -.016 0.0001 0 

17 29 17826.5 4.251 -.098 . 0044 0 

18 1300 40 146. 9 2. 167 0. 892 . 042 0 . 00086 

19 1300 37 335. 2 2. 525 . 874 . 008 0 . 00003 

20 1300 30 1245.7 3.095 . 845 f -.083 .0032 0 

21 1350 30 282.4 2.451 .658 0.329 -.083 . 0032 

22 1400 30 49.3 1.693 0 . 915 -.083 . 0032 

23 1400 22 339.0 2.530 0 .873 -.217 .0202 

24 1400 15 2287.2 3.359 j 0 .832 j -.384 | . 056 | 

k •> Pa)[A^ log t + II = PIT) • A • log t + pg) since AKON - 0 
r {o/Oo)° 1 _log(o/ao)' 

92 log <o/Oq) 93 4 ^ ^^^2 a In 10 

Figure 4. - Example of input matrix for MEGA16. Material is A-286 steel (ref, 1). Value of A =-0.05; 
AKON = 0; spline point =■ 1.56. Temperature stations are 1G5(]P, 1150P, 120QP, 1300° and 1400° F. 


H * P<T> + PCX) - G<SIGMA> 

33 ASTRGLOY 

0.000 -0.050 0.150 

0.121 0.107 0.079 


(a) Short printout. 


RESULTS FROM MEGIAB 

LDQH + A LDGH * PtT> f PtT 

G STATIONS 

I.IAO 1.&30 2.004 

P STATIONS 

1 400 . OOO 1 500 . OOO 1 600 . OOO 

B COEFFICIENTS 

2.373 - 4 . 327 -O . 706 

P STATION FUNCTION 

-3 . 204 - 1 . 502 O . OOO 

DATA USED IN CALCULATION 
LOG TIME 

OBS PRED DIFF 

1.10721 1. 10545 0.00176 

1 .770B5 1.89317 -0.12232 

2.24699 2.22124 0.02575 

2 . 60282 2 . 55466 0.04815 

2 . 76118 2 . 77860 -O .01742 

3 . 35790 3 . 27931 O . 07858 

3.60887 3.60017 0.00870 

1 . 48430 1 . 42398 O . 06032 

2.15290 2.10151 0.05139 

2 . 54568 2 . 57996 -O . 03428 

2 . 8S24B 2 . B04S4 O . 04794 

3 . 08930 3 . 1 4B03 -O . 05872 

3 . 34780 3 . 37789 -O . 03009 

3 . 64280 3 . 72627 -O . 08347 

1 .02119 0.98295 0.03824 

1 . 45939 1 . 53154 -O . 07214 

2. 16376 2.17941 -0.01565 

2.40312 2.44727 -0.04415 

2.72892 2.63887 0.09005 

2.94841 2.96600 -0.01759 

3.46235 3.39594 0.06642 

3.80147 3 . B5515 -O . 05367 

1 . 06070 1.20350 -O . 142GO 

1 . 64542 1 . 678 17 -O . 03275 

2.08243 2.01334 0.06909 

2.53491 2.44908 0.08583 

2 . 87315 2 . 75423 O . 11892 

3 . 24765 3 . 1 6826 O . 07940 

3.45312 3.59262 -0.13950 

0.7B533 0.81155 -0.02622 

1.69285 1.71094 -0.01809 

2.24055 2.16821 0.07234 

2 . 53237 2 . 55450 -O . 022 1 2 

33 ASTRQLOY ME 


> X G<SIGMA} 

33 ASTROLGY MEGA16 


1700.000 1800.000 

-46.084 

1.494 2.460 

THE VALUE OF 'A' - -0.150 

TEMP LOB STRESS G P 

SUM SO 

O . OOOOO 1400.0 2 . 00432 -1.71 770 - 3 . 20422 

0.01496 1400.0 1.93450 -0.65844 -3.20422 

0.01563 1400.0 1.90309 -0.21727 -3.20422 

O . 01795 1400 . O 1 . 86923 O . 23109 -3 . 20422 

O . 01825 1400 . O 1 . 84510 O . 53222 -3 . 20422 

O . 02443 1400 . O 1 . 78533 1 . 20554 - 3 . 20422 

0.02450 1400.0 1.74036 1.63700 -3.20422 

0.02814 1500.0 1.87506 0.15600 -1 .58200 

O . 03078 1500 . O 1 . 806 18 O . 98295 - 1 . 58200 

0.03196 1500.0 1.74819 1.56691 -1.58200 

0.03425 1500.0 1.71600 1.84101 -1.58200 

0.03770 1500.0 1.65321 2.26024 -1.58200 

0.03861 1500.0 1 .59106 2.54079 -1.58200 

0.04557 1500.0 1.49136 2.96600 -1.58200 

0.04704 1600.0 1.80618 0.98295 0.00000 

0.05224 1600.0 1.75205 1.53154 0.00000 

0.05249 1600.0 1.66745 2.17941 0.00000 

0.05444 1600 . O 1 . 61278 2 . 44727 O . OOOOO 

O . 06254 1600 . O 1 . 56820 2 . 63887 O . OOOOO 

0.06285 1600.0 1.49136 2.96600 0.00000 

O . 06726 1 600 . O 1 . 38917 3 . 39594 O . OOOOO 

O . 070 15 1 600 . O 1 . 27875 3 . 85515 O . OOOOO 

O . 09054 1700 . O 1 .61278 2 . 44727 1 . 49449 

0.09161 1700.0 1.52504 2.82306 1.49449 

0.09638 1700.0 1.46240 3.08840 1.49449 

O. 10375 1700.0 1.38021 3.43337 1.49449 

0.11789 1700.0 1.32222 3.67494 1.49449 

0*12419 1700.0 1.24304 4.00272 1.49449 

0.14365 1700.0 1.16137 4.33868 1.49449 

O . 1 4434 1800 . O 1 . 46982 3 . 05707 2 . 46033 

0.14467 1800.0 1.31175 3.71840 2.46033 

0.14990 1800.0 1.23045 4.05464 2.46033 

0.15039 1800.0 1.16137 4.33068 2.46033 

:GA16 

THE VALUE OF 'A' ■ -0.150 


S. D. OF RES X 0.079 DOF - 24.0 SUM OF RES - 0.012 AV ABS DEV - 0.057 


(b) Long printout. 

Figure 5. - Example of output from MEGA 16. 
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ASTROLOY MEG16B AK0N = 3-8r 

A • -B . ) 5 



LOG TIME. HOURS 
(a) Plot of input data and predicted isothermals. 


ASTRaLO> 


MEG16B 


AKOK'= 30 

ft = -0.1b 



(d) Plot of temperature function. 


33 A5TRGLCY KEG16B AKON'-30 

ft = -e.ib 



; i I ! I I < 

.p I.? l.*» !.4 l‘B 

LOG STRESS. <SI 
(b) Plot of stress function versus log stress. 



ASTRCLCV KEGloB AKO\-30 

ft -0.1b 


33 ASTRG^OY i^lGISB AkGR=3C 

A = -B.lb 


2 .b - 


-B. 32 ZE 01 =• Pj^ at lowest temperature 
- 0 .I 58 E 01 = ?2 at second temperature 
0.a“-bE 01 ' P4 at highest temperature 


2.hL 0 .P 37 E 01 = Yi ' function at spline point 



tc) Plot of stress function versus stress on linear scale. 


(f) Plot of the "master curve*. 


Figure 6. - Sample plots. 



STRESS* KSI _ OBS 



33 i^STROLOY 


MEG16B AK0N = 3£f 

A • 0.00 




e' i - J J i . _ 1 I I I I 1 

B-5 0 .’> l.J 1.7 S.I 2,5 2.9 3.3 3.7 9.1 9,5 

LOG TIME, hours 

(h) Plot of raw data for preassessment using linear stress scale. 

,, 33 ASTROLOY ALL MLG]6B 


o 

ift 


.* L 4 I 4- 4 I I 1 I 

-.iM -.15 -.!• -.«s «.«• «.es a. IS a. IS 0 . 2 a 

A VALUE 

(i) The A scan: observed standard deviation versus selected A value. 

Figure 6. - Concluded. 


rapidly to about 5 seconds per data set, if many sets 
are run consecutively after the program is loaded. 
When run on this virtual memory computer MEGAI6 
uses approximately 2800 K of eight bit bytes. 


Discussion 

The question of accuracy, somewhat exemplified 
by the observed versus predicted plot in each run of 
MEGAI6, has been discussed (refs. 1 to 3), but not 
really answered. A point not always considered is 
that the accuracy of a prediction cannot be much 
better than the precision of the data used to make 
that prediction. Of course the scatter of a given set of 
creep-life data depends on many unquantified 
variables in the general categories of test procedure, 
metallurgical differences, environment, etc. 
Estimates of the extent or importance of this scatter 
are not often presented with creep-life data. Few 
researchers are inclined to run the expensive replicate 
tests necessary to given meaningful estimates of the 
experimental error for various materials and 
laboratory conditions. To provide a rough estimate 
of scatter to be found in creep-life experiments, the 
data analyzed in the development of mega 1 6 (ref. 4) 
were also used to provide the measures of average 
deviations in table I. Shown are some statistics for 
the temperature and stress combinations which had 
three or more repeat tests. Using the usual 
assumption of a lognormal distribution of time to 
failure, the last column of this table shows a fairly 
wide variation in the scatter about the mean log life 
for various combinations of material temperature 
and stress. For some materials standard deviations of 
0. 1 appear to be typical, but for others 0.2 or over are 
to be expected. Therefore, it is reasonable to assume 
that the best accuracy to be achieved from a 
predictive technique will also be somewhere in the 
range of 0. 1 of a decade of the logarithm of the time 
to rupture, but there will be cases of higher deviation 
for certain materials and testing conditions. 

With the stress function of MEGAI6, there is the 
need for choosing an intermediate value of the 
logarithm of stress to be used as the spline point. 
Normally, a value somewhere near the midrange of 
stress will suffice, but for some data selection of 
another log stress value will change the fit and give 
different extrapolations. As an example, some 
additional calculations were made using the type 316 
stainless steel set of data used in reference 4. In that 
report the location of the spline point for this 
material had been, for other reasons, selected at a 
point somewhat distant from the midrange. Figures 7 
and 8 show the isothermals and master curves for 
some of the selected combinations of A and spline 
point locations. Table II gives the resulting standard 
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TABLE I.— STATISTICS OF REPEAT MEASUREMENTS 


TEMP 

STRESS 

NO* OF 

LIFE 


LOG 

F 

KSI 

REPEATS 

HEAN 


C.V. 

MEAN 




(a) 


(b) 

(a) 




64 IlOO 

-O ALUMINUM 

482* 

6*0 

3 

299*67 


41*63 

2.445 

482. 

5*5 

3 

1274.67 


58*67 

3*055 

662. 

4*0 

3 

49.17 


7.22 

1.691 

662* 

3.5 

4 

203*75 


14.28 

2*306 

752* 

3*5 

6 

13.08 


37.86 

1.097 

752* 

3.0 

7 

62*71 


15.99 

1.793 

932. 

1.5 

3 

224*33 


26*17 

2*340 




75 5454- 

-Q ALUMINUM 

662* 

14.0 

3 

81.67 


26*67 

1.902 

662. 

11*0 

5 

436*00 


14*32 

2*636 

752. 

9*0 

7 

165.71 


13.50 

2.216 

842* 

7*0 

4 

96.50 


9.48 

1*983 

932. 

5*0 

3 

113.00 


13*74 

2*050 




95 INCO 

625 


1200* 

70*0 

6 

46*57 


83*72 

1*501 

1200. 

55*0 

3 

769*53 


69*57 

2*781 

1300* 

50*0 

4 

46*67 


80.86 

1.549 

1300. 

45*0 

7 

246.94 


54*79 

2*307 

1400. 

30.0 

5 

267*52 


69*10 

2*326 

1500. 

20.0 

5 

43.18 


32*88 

1*617 

1500. 

17*5 

5 

87*66 


22.75 

1.933 

1500. 

15.0 

9 

506*80 


73*06 

2*601 

1500. 

12*5 

4 

1 722 . 42 


53.93 

3*171 

1500* 

10.0 

5 

7117*70 


34*85 

3*834 

1600* 

12*0 

3 

44.00 


19*38 

1*638 




105 U- 

-500 


1200. 

140*0 

6 

7.53 


35*93 

0*846 

1200* 

100*0 

8 

1092*57 


57.37 

2.981 

1500. 

60*0 

8 

17*99 


31.58 

1*235 

1500* 

30*0 

6 

1060*93 


50*99 

2*959 

1800* 

9*0 

3 

83.57 


59*33 

1*849 




104 L- 

605 


1200* 

60.0 

8 

7*56 


38*03 

0*851 

1200* 

40*0 

6 

1129*18 


30.59 

3*035 

1500* 

30*0 

8 

14*45 


28*42 

1.147 

1500* 

17.0 

7 

1618*47 


19.68 

3.201 




lOl 6061- 

-T6 ALUMINUM 

662* 

26*0 

3 

173*00 


18*05 

2*234 

752. 

26*0 

3 

28.37 

133.22 

1 . 163 

752* 

21.0 

6 

76*67 


20.27 

1.878 

842. 

13.0 

4 

184*38 


30*40 

2*250 

842. 

11*0 

3 

751.33 


22* lO 

2*869 

932* 

13*0 

4 

19.30 


59.09 

1.225 

932* 

11.0 

3 

74*33 


12*50 

1.869 

1112* 

4*0 

3 

133.33 


7.09 

2.124 


^MEAN= — 
n 


^C.V. 


VvAR 

MEAN 


xlOO 


S,D. = V VAR where VAR = 


ejc2-(i:x)V/? 


LIFE 

S*D. 

(C) 


0«212 

0.258 

0.031 

0.061 

0.136 

0.070 

0*122 

0.127 


0.112 

0*063 

0*059 

0*041 

0*060 

0*067 


0.476 

0*409 

0*382 

0*336 

0*350 

0*139 

0*107 

0.321 

0*298 

0*139 

0*082 

0*276 


0*192 

0*239 

0*141 

0*288 

0*337 

0*239 


0*166 

0*140 

0*111 

0*092 

0.127 


0*076 

0*602 

0.077 

0.134 

0*092 

0.267 

0.056 

0*030 

0.167 


n-\ 






Figure 8. - Results from MEGA16 showing effect of different A values and 
spline points on shape of master curves. Material, 316 stainless steel. 





Figure 7. - Results from MEGA16 showing effect of different A values and spline points on shape 
of predicted isothermals. Material. 316 stainless steel. 


TABLE II.— EFFECT OF 

LOCATION OF SPLINE 

POINT IN THE STRESS 

FUNCTION ON THE 

STANDARD DEVIATION 

OF PREDICTIONS 

[An example using type 316 
stainless steel with 38 data 
and a stress range of 7 to 
30 ksi. The program is 
MEGA 16 with AKON = 30.] 


Location 
of spline 
point 

A value 

-0.05 

0 

0.05 

Standard deviation 

1.0 

0.130 

0.125 

0.125 

1.1 

0.123 

.118 

.116 

1.2 

0.109 

.103 

.102 

1.3 

0.082 

.076 

.076 

1.38 

0.082 

.080 

.082 

1.4 

0.085 

.083 

.085 
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deviations for all the combinations. From the table, 
it appears that a minimum standard deviation exists 
near 0.076, which is a lower value than that reported 
in figure 39 of reference 4. Thus, choice of spline 
point location might be a means of fine tuning the 
predictions of mega 1 6 , but from figure 7(a) it is 
obvious that extreme values of the spine point may 
give strange predicted isothermals. 


Summary of Results 

MEGA16 is the outcome of a number of years of 
study of various forms of the so-called minimum 
commitment method of analysis of stress-rupture 
data. Basically, it is a Fortran IV program which fits 
a specific, yet flexible, model equation to a given 
family of time-to-failure data. The model includes a 
discrete temperature function, a two-part stress 
function, and an adjustable constant. A, which 
affects the extrapolation beyond the range of data. 
MEGA16 can handle up to 200 data points per set, 
including the weighting of certain points, if so 
desired. For each datum it sets up an equation 
according to the assumed model and some optional 
locations of the centers of the temperature and stress 
functions. After it solves the resulting set of 
equations to determine the coefficients of the terms 


of the model equation, it predicts the rupture life for 
each input combination of temperature and stress 
and predicts the logarithm of the time of evenly 
spaced values of the logarithm of stress. 

MEGA16 produces many different types of plots to 
depict the response of a particular material. These 
include (1) the basic predicted time-to-rupture versus 
stress for specific temperatures, (2) the temperature 
and stress functions, (3) the time and temperature 
function versus the logarithm of stress (the master 
curve), (4) the observed versus predicted log time to 
failure values, (5) residual plots of the logarithm of 
time, stress, and temperature, and (6) standard 
deviation as a function of the parameter A. MEGA16 
calculates, as an estimate of accuracy, the standard 
deviation of the observed minus the predicted 
rupture life. 

To be analyzed by megai6, the data set 
representing a material should include test results 
from at least three different temperatures, each with 
several levels of stress, and the resulting spread of life 
values. That set can be run alone or with many other 
sets of data from other materials to obtain complete 
analyses, in easy-to -understand graphical form. 

Lewis Research Center 

National Aeronautics and Space Administration 
Cleveland, Ohio, August 28, 1980 
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Appendix— Listing of MEGA16 and its Subroutines 


THIS IS M E G A I 6 THE MANSON-ENSIGN GENERALIZED ANALYSIS 

OF STRESS RUPTURE DATA INCORPORATING THE IDEAS OF HINIMUM 

COMMITTMENT TO FORM OF EQUATION, THE STATION FUNCTION APPROACH, 

AND A MODEL EQUATION OF THE FORM 

LOGH + (1 + A LOGH)P(T) □ G(S) 

WHERE H = TIME 

P = STATION FUNCTION FOR TEMPERATURE 
T = TEMP 

G = STRESS FUNCTION, D + E LOG S + F SKXALPHA 
ALPHA = +1 OR -1 
S = STRESS / REFERENCE STRESS 
A = A CONSTANT USUALLY BETWEEN -.2 AND +.2 
OR A = A(1-ACON(T-TMID/TMID+A60))>«K2 


ASSUMPTIONS 

FORM OF THE MODEL EQUATION 

STRESS AND TEMPERATURE ARE THE ONLY 
INDEPENDENT VARIABLES 


RELATIONSHIPS AMONG STRESS RUPTURE LIFE, 
STRESS, AND TEMP PRESENT IN THE SHORT 
TIME OBSERVATIONS ALSO EXIST IN 
THE LONG TIME BEHAVIOR 


CALL STRUCTURE 


MEG16B 


TU5TR 


STACN 


AMATF 

DMFSS 


DML5S 

TOPLT 

CVFITF 


DTRIA 

DSOLVE 

FUNC 
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C THE SUBROUTINES CALLED BY THIS PROGRAM (MEG16B) ARE 

C 

C 

C NAME SOURCE* PURPOSE 

C 

C TUSTR STU5TR CONVERTS LOG STRESS VALUE TO STATION FUNCTION FORM 

C STACN SSTACN CONVERTS TEMP VALUE TO STATION FUNCTION FORM 

C AMATF SAMATF SETS UP THE MATRIX OF STATION FUNCTION EQUATIONS ( UMATRX) 

C TOPLT STOPLT PREPARES STATION FUNCTION VALUES (XF) FOR PLOTS AND CURVE FITS 
C 
C 
C 

C OTHER SUBROUTINES CALLED ARE 

C 


C 

c 

IN 

NAME 

SOURCE. 

PURPOSE 




L 

c 

AMATF 

DMFSS 

SDMF 

PREPARES MATRIX, 

FINDS 

RANK, 

AND LINEARLY IND. ROWS 

C 

AMATF 

DMLSS 

SDML 

INVERTS UMATRX TO 

FIND 

LEAST 

SQUARES SOLUTION (XF) 

c 

TOPLT 

CVFITF 

CRVFF 

POLYNOMIAL CURVE 

FIT — 

CROUSE-JORDAN REDUCTION 

c 

CVFITF 

DTRIA 

DCTRIA 





c 

CVFITF 

DSOLVE 

DCSOLV 





c 

CVFITF 

FUNC 

FUNCT 






C 

C ADDITIONAL SUBROUTINES CALLED ARE FROM THE IBM PLOTTING PACKAGE 

C 

C THESE INCLUDE GPLOT, XAXIS. YAXIS, CHARS, NUMBER, BEGID, TITLE, SCLBAK, ETC. 

C 
C 
C 

C NOMENCLATURE INPUT VALUES 

C 

C TITEL IS THE TITLE FOR A PARTICULAR MATERIAL 
C NF DENOTES THE F COLUMN TO BE OMITTED IN FINAL MATRIX 

C NP DENOTES THE P COLUMN TO BE OMITTED IN FINAL MATRIX 

C NG DENOTES THE G COLUMN TO BE OMITTED IN FINAL MATRIX 

C MK DENOTES THE COLUMN TO BE USED AS THE DEPENDENT VARIABLE 

C NL IS THE NUMBER OF STATIONS FOR TIME— NOT USED IN THIS VERSION-SET = NT 
C NT IS THE NUMBER OF STATIONS FOR TEMP MAXIMUM = 32 

C NS IS THE NUMBER OF STATIONS FOR STRESS MAXIMUM = 16 

C NQ IS THE NUMBER OF STRESS PARAMETERS BEING DETERMINED 
C NCOL IS THE NUMBER OF COLUMNS IN UMATRIX = NUMBER OF BOTH 

C DEPENDENT PLUS INDEPENDENT VARIABLES MAXIMUM = ^8 

C NEV IS THE NUMBER OF EVENLY SPACED LOG TIME VALUES AT 
C WHICH PREDICTIONS ARE MADE 

C NISO IS THE NUMBER OF ISOTHERMALS TO BE PREDICTED 
C WHICH ARE NOT ALSO STATION FUNCTIONS 

C NXAV IS NUMBER OF EXTRA A VALUES READ IN (MAX=1 FOR NOW) 

C CISOT(K) ARE THE TEMP VALUES OF THE ISOTHERMALS 

C ALOGT ARE THE LOG TIME STATIONS (NOT USED IN THIS VERSION) 

C TSTA ARE THE TEMP STATIONS 

C S5TA ARE THE LOG STRESS STATIONS 

C NODU IS THE NUMBER OF INPUT DATA POINTS — DATA USED 

C NODUW IS THE NUMBER OF DATA USED AFTER WEIGHTING 

C NODP IS THE NUMBER OF LONG TIME DATA TO BE PREDICTED (NOT IN THIS VERSION) 

C NOEXT IS THE NUMBER OF DATA ADDED BY THE WEIGHTING 

C CUTOF IS THE CUTOFF LIFE IN HOURS 

C AEXT IS THE EXTRA A VALUE READ IN WITH DATA 

C OPT IS THE PLOT OUTPUT OPTION; F (OR 0) MEANS PLOT, T (1) MEANS NO PLOT 

C CONSTR ARE THE VALUES OF CONSTANT STRESS FOR THE 

C ISOSTRESS PLOTS 

C TMP(I) ARE THE TEMP VALUES OF THE DATA IN DEGREES FAHREHEIT 

C SIG(I) ARE THE STRESS VALUES OF THE DATA IN KSI 

C ALS(I) ARE THE LOG STRESS VALUES 

C AT(I) ARE THE TIME TO FAILURE VALUES OF THE DATA IN HOURS 

C ALT(I) ARE THE LOG TIME VALUES 

C WF IS THE NUMBER OF TIMES A POINT IS TO BE WEIGHTED 
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NOMENCLATURE xxxxx CALCULATED VALUES 

PREFIX W USUALLY MEANS WEIGHTED VALUE 

PREFIX Q USUALLY MEANS SINGLE PRECISION FOR PLOTTING 

XF IS THE SOLUTION VECTOR FROM THE INVERSION OF THE SET OF EQUATIONS 

SFCN(J) ARE THE STRESS STATION FUNCTION VALUES— THE STRESS FUNCTION 
TFCN(J) ARE THE TEMP STATION FUNCTION VALUES— THE TEMP FUNCTION 

STDEV IS THE STANDARD DEVIATION OF THE RESIDUALS (THE DATA USED) 

STDL IS THE STANDARD ERROR OF THE PREDICTED VALUES = RMS VALUE 

AVARSL IS THE AVERAGE ABSOLUTE DEVIATION OF THE PREDICTED VALUES 


NOTA BENE MEG16B 

MAXIMUM NUMBER OF INPUT DATA POINTS THIS VERSION CAN TAKE IS 200 
MAXIMUM SIZE OF MATRIX TO BE INVERTED IS AS X AS AMATF,DMFSS 

DOUBLE PRECISION SSTA(16), C0EFG(16), C0EFGG(16) 

DOUBLE PRECISION SFCN(16) 

DOUBLE PRECISION AL0GTC32), COEFPC32), TSTA(32), TMPTC32) 

DOUBLE PRECISION DAZERO(AS), DTFIL(6A) 

DOUBLE PRECISION ALS(200), ALTC200), TRIP(200) 

DOUBLE PRECISION AT(200), PRDD(200), PRODC200), SIG(200) 

DOUBLE PRECISION TMP(200), WALS(200), WALT(200), WTMP(200) 

DOUBLE PRECISION UMATRXC 200 , AS ) , VALUA(S), VALUM(200 ) 

DOUBLE PRECISION XF, XPLTT, TFCN, XPLTRT 

DOUBLE PRECISION XPLTG, YPLTG, SCOEF, AEXT 

DOUBLE PRECISION TVALP, YVALU, YVAL, YVALG 

DOUBLE PRECISION SDEL, SVAL, SVLG 

DOUBLE PRECISION Cl , C2 , C3 , CA , C5 , UNKN( 1 1 ) 

DIMENSION QLABX3(1), QLABX7(1), QLABY2(1), QLABY3(1) 

DIMENSION QLABY7(1), QPX(2), QPX2(2), QPXL(2) 

DIMENSION QPXL2(2), QPXR(2), QPXR2(2), QPYC2) 

DIMENSION QPY2(2), QPYL(2), QPYL2(2), QPYR(2) 

DIMENSION QPYR2(2), QLABX2(3), QLABXD(A) 

DIMENSION QLABXKA), QLABYD(A), QLABYKA) 

DIMENSION QABKEY(5), QADKEY(5), QAVKEY(5), QEXKEY(S) 

DIMENSION QRDKEY(5), QSDKEY(S), QLABX5C7), QLABXAC7) 

DIMENSION QLABXC(S), QLABYC(S), PSDEV(9), PSTDEVC9) 

DIMENSION PSTDL(9), C0NSTR(I6) 

DIMENSION TITEL(16), CISOT(32), PIS0C32) 

DIMENSION ALOGP(AS), ALSPE(AS), ALTP(AS) 

DIMENSION ALTPE(AS), ALTTC3A), QLABXA(2), QLABYA(5) 

DIMENSION AZERO(AS), BNONZ(AS) 

DIMENSION ALSTRLC6A), PFILC6A), TFIL(6A) 

DIMENSION TRFILC6A) 

DIMENSION RESL(96), ABRES(200), GCALCC200), HOUR(200) 

DIMENSION HOURS(200), QALS(200), QALT(200), QXPL(200) 

DIMENSION QYPL(200), RES(200), RESD(200) 

DIMENSION RSDEV(200). UF(200), PEST(200) 

DIMENSION AMATRXC AS, 9),BMATRX(A8,9),QVALUA(8) 

DIMENSION QANGLS(16), QLABXS(3), PLRESD(200), QLABYA(5) 

DIMENSION QTMPL(200), QSIGL(200), QTMP(200) 

DIMENSION ZDEV(200), PSDZD(9), XPLTRS(6A), GPLT(6A), ALGTC(6A) 

DIMENSION GPLAT(200), STRLC6A), GLOWC6A), GHIGH(6A) 

DIMENSION IVGR(7) ,GRAXV(10) ,GRAYV(10) ,KSYM(30) 

DIMENSION GRLEG5(2) ,GRLEG6(A) ,XMIS(2) ,YMIS(2) 

LOGICAUa XAXV,YAXV 
DATA XAXV /.TRUE./ 

DATA YAXV /.FALSE./ 

INTEGER)<2 NUMPT/20/ 

LOGICAU«l OPT(17) 
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c 


COMMON 

COMMON 

COMMON 

COMMON 


/ACOMN/ QXPLTTC32),QXPLRTC32),QYPLTT(52),QXPLTGCI6),QYPLTG(16) 
/BCOMN/ XPLTTC32) ,TFCN(32) ,XPLTRT(32) ,XPLTG(16),YPLTG(16) 
/CCOMN/ XF(48),NP 
/DCOMN/ D,E,F,G,H,AA 


PRINT OUT TODAYS DATE ON EACH PLOT 
INTEGER^2 NDATE 

DIMENSION DATES(<t) 

DATA NDATE /16/ 

DATA DATE5(1),DATES(2) / ’DATE*,* — » / 
DATA DATESC3) ,DATES(4) / ’MO/DS’D/YY» / 
CALL TIME (NDATE, DATES) 


WRITE (7,516) 

(DATES(I), 1=1,4) 



DATA 

QLABXl 

/ 

» LOG 

• , ’ TIM’ , ’E, H’ 

, ’OURS’/ 


DATA 

QLABYl 

/ 

’ LOG 

’,’ STR’,’ESS,’ 

,’KSI ’/ 


DATA 

QLABY2 

/ 

'G(S) 


/ 


DATA 

QLABX2 

/ 

’LOG 

’,’STRE’,’SS ’ 

/ 


DATA 

QLABX3 

/ 

’TEMP 


/ 


DATA 

QLA3Y3 

/ 

’P(T) 


/ 


DATA 

QLABY4 

/ 

’RESI 

’,’DUAL’,’S O’ 

, ’F L’ , ’OG T’ 


DATA 

QLABX5 

/ 

’LOGH 

’ , ’ + A’ , ’>^LOG’ 

,’H?(P(’,’T) +’,’ P(T’,’) 

ff 

DATA 

QLABX7 

/ 

’PRED 


/ 


DATA 

QLABY7 

/ 

’ OBS 


/ 


DATA 

QLABX8 

/ 

’STRE 

’,’SS, ’,’ KSI’ 

/ 


DATA 

QLABX4 

/ 

1 

’ , ’1000’ , ’ / T’ 

, ’EMP ’ , ’+ 46’ , ’0 ’ , ’ 

’/ 

DATA 

QLABXA 

/ 

’A V’ 

, ’ALUE’ / 



DATA 

QLABYA 

/ 

’S.D. ’ 

,’ OF ’,’PRED’, 

’ICTI’,’ONS ’ / 


DATA 

QLABXD 

/ 

’LOG 

’ , ’TIME’ , ’ ,A =' 

, ’ 0 ’ / 


DATA 

QLABYD 

/ 

’LOG 

’ , ’TIME’ , ’ ,A =’ 

, ’ X ’ / 


DATA 

QLABXC 

/ 

’LOG 

’ , ’T, A’ , ’=0 ’ 

, LO’ , ’G CU’ , ’TOFF’ , ’ 


DATA 

QLABYC 

/ 

’LOG 

’,’T, A’,’=X ’ 

LO’,’G T,’,’ A=0’,’ 

' , ’ 


DATA KSYM /7 0 , 66 , 62 , 65 , 68 , 56 , 6 7 , 6 1 , 72 , 92 , 78 , 6 9 , 58 , 6 3 , 7 1 , 

1 70,66,62,65, 92,63,67,61,72,78,69,71,56,1^3,181 / 

DATA NEV /3A/, ALTT / 1.0,1.125,1.25,1.375,1.5, - 

11.625.1.75.1.875.2.0. 2.125.2.25.2.375.2.5.2.625.2.75.2.875.3.0, 

23.125.3.25.3.375.3.5.3.625.3.75.3.875.4.0. 4.125.4.25.4.375.4.5, 
34.625,4.69897,4.75,4.875,5.0/ 


NOTA BENE 


C 


C 


THE FIRST A SHOULD BE ZERO TO INSURE THAT ALL PLOTS HAVE MEANING 
SET LAK TO AT LEAST 1 LESS THAN DIMENSIONS OF VALUA 


DATA 

LAK /5/ 

, VALUA /0.0,-.15 

,-0.1,- 

0.05, .15/ 


DATA 

LAK /2/ 

, VALUA /O. 00,-0. 

05/ 



DATA 

GRLEG5 

/ ’ A ' ’ 

/ 



DATA 

QPX /O 

.5,5.3/ 




DATA 

QPY /O 

.5,5.3/ 




DATA 

QSDKEY 

/ ’SD O’ , ’F RE’ , 

’S =’, 

» f f 

, 

» / 

DATA 

QEXKEY 

/ ’RMS ’ , ’OF P’ , 

’RED=’ , 

f f f 

, 

’/ 

DATA 

QAVKEY 

/ ’AV D’ , ’EV P’ , 

’RED=’, 

^ * 

'/ 

DATA 

QABKEY 

/ ’ AB D’ , ’EV P’ , 

’RED=’ , 

1 V T 

, 

’/ 

DATA 

QADKEY 

/ ’AV A’ , ’BS D’ , 

’EV =’ , 

1 f » 

, 

’ / 

DATA 

QRDKEY 

/ ’SD R’ , ’EDUC’ , 

’ED =’, 

f r r 

, 

’ / 


DATA QPXL 
DATA QPYL 
DATA QPXR 
DATA QPYR 
DATA QPX2 
DATA QPY2 
DATA QPXL2 
DATA QPYL2 
DATA QPXR2 
DATA QPYR2 


/O. 5, 4. 699/ 
/O. 801, 5.0/ 

/O .801,5.301/ 
/0.5,5.0/ 

/2. 0,5.0/ 

/2. 0,5.0/ 

/2. 0,4.699/ 
/2. 301, 5. 0/ 
/2. 301, 5.0/ 
/2. 0,4.699/ 
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ir 


INITIALIZE— -BEGIN WITH A NEW DATA SET 

100 DO 105 1=1,200 
TMPC I) =0 . 0 
SIG(I)=0.0 
AT(I)=0.0 

105 WF(I)=0.0 

DO 106 1=1,9 
PSTDEV(I)=0.0 
PSDEV(I)=0.0 
P5TDL(I)=0.0 

106 P5DZD(I)=0.0 
5TDEV = 0.0 

DO 107 J=l,200 
ALTCJ) = 0.0 
ALSCJ) = 0.0 

107 CONTINUE 

READ (5,^10) TITEL 

INPUT-STATION VALUES AND CONSTANTS 

READ (5,A02) NF,NP,NG,MK 

READ (5,402) NL,NT,NS,NISO,NXAV 

READ (5,400) ( CISOT ( K ) , K=1 , NISO ) 

READ (5,400) ( A L OGT ( I ) , I =1 , NL ) 

READ (5,400 ) ( TS T A ( I ) , I = 1 , NT ) 

READ (5,400) ( S5TA ( I ) , I = 1 , NS ) 

READ (5,404) NODU , NODT , CUTOF , AEXT , ( OPT ( I ) , I = 1 , 17 ) 

NCS=8 

READ (5,400) ( CONSTR ( K ) , K=1 , NCS ) , ALPHL , ALPHH 
WRITE (6,515) 

WRITE (7,515) 

WRITE (9,515) 

WRITE (6,455) TITEL 
WRITE (7,455) TITEL 
WRITE (7,513) 

WRITE (9,455) TITEL 

WRITE (9,463) ( SST A ( I ) , I =1 , NS ) 

WRITE (9,467) ( TSTA ( I ) , I =1 , NT ) 

WRITE (7,477) NF , NP , NG , MK , NL , NT , NS , NISO , NXAV 
KTEST=1 

CUTLOG ° ALOGIO(CUTOF) 

CUTHAF ° (CUTLOG+5.0)/2.0 

DO 111 1=1, NS 

QANGLS(I) □ lO.C^^^SSTAd) 

111 CONTINUE 

IF (NXAV.NE.l) GO TO 109 
LAK = LAK+MXAV 
VALUA(LAK) = AEXT 
109 NOEXT = 0 

INPUT--OBSERVATIONS 

DO 115 1=1, NODT 

READ (5,408,END=399,ERR=113) TMP ( I) , SIG( I ) , AT ( I) , WF( I ) 
ALT(I)=DLOG10CAT(I) ) 

ALS(I)=DLOG10(SIG(D) 

GO TO 114 

113 CONTINUE 

WEIGHT DATA 

114 IF (NOEXT.GT.O) GO TO 117 
KN = I 

GO TO 118 
117 KM = KN+1 
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lis WTMPCKN) = TMP(I) 

QTMP(KN) = TMPC I) 

WALSCKN) = AL5CI) 

QAL5CKN) = ALSCn 
WALTCKN) = ALT(I) 

QALT(KN) = ALT(I) 

IF (WFCDAE.l.O) GO TO 116 
KUF = WFCI)“1.0 
DO 112 KK= 1,KWF 
KH = KN+1 
WTMPCKN) = TMPCI) 

WALSCKN) = ALSCI) 

QAL5CKN) = ALSCI) 

WALTCKN) = ALTCI) 

NOEXT = NOEXT+1 
112 CONTINUE 

116 IF CI.EQ.l) GO TO 115 

IF CWTMPCI) .EQ.WTMPCI-D) GO TO 115 
KTEST=KTE5T+1 
115 CONTINUE 

NODUW □ NODU+NOEXT 
NODTW □ NODT+NOEXT 

CHOOSE AN A VALUE 

121 DO 385 LA=1,LAK 

INITIALIZE ALL THE VECTORS IN COMMON 
DO 119 K=l,^8 
XFCK) = 0.0 
XPLTTCK) = 0.0 
TFCNCK) = 0.0 
XPLTGCK) = 0.0 
YPLTG(K) = 0.0 
QXPLTTCK) □ 0.0 
QYPLTTCK) = 0.0 
QXPLTGCK) = 0.0 
QYPLTGCK) □ 0.0 
QXPLRTCK) = 0.0 
XPLTRTCK) = 0.0 

119 CONTINUE 

DO 123 1=1,16 
SFCNCI)=0.0 
123 CONTINUE 

DO 120 1 = 1, A8 
DO 120 J=l,200 

120 UMATRXCJ,I)=0.0 

BEGIN PLOT OF ISOTHERMALS 

CALL BEGID(l) 

IVGRCl) = 7 
IVGRC2) = 1 
IVGRC3) = 3 
IVGRCA) = 70 
IVGRC5) □ 1 
IVGRC6) = 15 
IVGRC7) = 0 
GRAXVCl) = 10. 

GRAXVC2) = -1. 

GRAXVC3) = 0. 

GRAXVCL) = 1. 

GRAXVC5) □ 5. 

GRAXVC6) = 8. 

GRAXVC7) = 2. 

GRAXVC8) = -1. 

GRAXVC9) = 0. 

GRAXVCIO) = 12. 

GRAYVCl) = 10. 
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GRAYV(2) □ -1. 

GRAYVC3) = 90. 

GRAYV(A) = .0 
GRAYVC5) = 2.5 
GRAYVC6) = 5. 

GRAYVC7) = 2. 

GRAYV(S) □ “1. 

GRAYVC9) = 0. 

GRAYV(IO) = 12. 

CALL XAXIS ( .9, .6,GRAXV) 

CALL YAXIS ( .9, .6,GRAYV) 

IA = 1 

DO 1^0 K=1,KTEST 
IK = 0 

DO 130 I=IA,N0DTW 
IK=IK+1 

IF (I.EQ.IA) GO TO 125 
IF (WTMP(I) .NE.WTMPCI-D) GO TO 135 
125 QXPL(IK)=UALT(I) 

QYPL(IK)=WALS(I) 

130 CONTINUE 
135 IA=I 

IPLT^IK’l 

IF (K.EQ.KTEST) IPLT=IPLT+1 

PLOT THE DATA FOR THE K TH TEMPERATURE 

IVGRC2) = IPLT 
IVGRC3) = 3 
IVGR(A) = KSYM(K) 

CALL GPLOT ( QXP L , QYP L , I VGR ) 
l<+0 CONTINUE 

IIRIIE (6,518) (IVGR(I),I = 1,7) 

WRITE (6,501) (GRAXV(I),I=1,10) 

WRITE (6,501) (GRAYVC I ) , 1 = 1 , 10 ) 

WRITE (6,519) ((QXPL(I),I=1,6),(QYPL(I),I=1,6)) 
CALL TITLE ( 1 , 6 A , 2 5 , T I T EL ) 

CALL CtOARS ( 1 6 , DA T ES , 0 . , 1 . 7 , 9 . , 1 2 ) 

CALL TITLE ( 3 , 1 6 , 2 0 , Q L ABYl ) 

CALL CHARS ( 8 , GR L EG5 , 0 . , 8 . 0 , 9 . 5 , 1 5 ) 

CALL NUMBER ( A , V A L U A ( L A ) , 8 , 2 , GR L EG6 ) 

CALL CHARS ( 8 , GR L EG6 , 0 . , 8 . 5 , 9 . 5 , 1 5 ) 

CALL TITLE ( ^ , 16 , 20 , QLABXl ) 

FORM SET OF SIMULTANEOUS EQUATIONS (UMATRX) 


NO = A 

AKON = 0 . 0 

lAl NCOL=NQ+l+NT 

IF (NP.GT.O) NCOL=NCOL“l 
NOTA BENE MEG16B 

NUMBER OF UNKNOWNS IS ^ + NT 1 NUMBER OF COLS IN UMATRX IS ^ + NT 

THIS USES S □ STRESS / STRESS MID 
ALSO ALPHH = EXPONENT ON S HIGH STRESS 

ALPHL = EXPONENT ON S LOW STRESS 

IF (ALPHH. NE. 0.0) GO TO 1^3 
ALPHH = -1.0 

ALPHL = 1.0 

1^3 DO 210 J=l,NODUW 
142 SCOEF=WALT( J) 

YVALG = WALS(J) 

TRIP( J)=( (WTMP( J)-TSTA(NP) )/(TSTA(NP)+460 . ) )^^2 
VALUrU J ) aVALUA( LA)X( 1 . 0 - AKON ^ TR I P ( J ) ) 

161 CALL TUSTR ( SSTA , YVALG, NQ , ALPHL , ALPHH , COEFG) 

C WRITE (7,450) ( COEFGC K ) , K=1 , NQ ) 
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162 KNT=0 

DO 170 K=1,NT 

IF (WTMP(J) .NE.TSTA(K)) GO TO 165 

163 PRODC J)=VALUMC J)^WALT(J)+1.0 
TMPT(K)=PROD( J) 

GO TO 170 
165 TMPTCK)=0.0 
KNT=KNT+1 

170 CONTINUE 

171 IF (KNT.NE.NT) GO TO 190 
C 

TVALP □ WTMPCJ) 

186 CALL STACN ( TST A , TVALP , NT , TMPT , N ) 

PRDDCJ) -VALUMCJ)^WALT(J)+1 . 0 

187 DO 193 K=1,NT 

192 TMPT(K) = TMPT(K)^PRDD( J) 

193 CONTINUE 

C WRITE (7, ^^5) (TMPT(K) ,K = 1.NT) 

190 CONTINUE 

194 N=1 

UMATRXCJ,M)=SCOEF 
DO 200 1=1, NT 
IF (I.EQ.NP) GO TO 200 
M=M+1 

UMATRXC J,M)=TMPT(I) 

200 CONTINUE 

DO 205 1=1, NQ 
N=M+1 

UriATRXC J,M)=COEFG(I) 

205 CONTINUE 
210 CONTINUE 

IEND=NODUW 

225 K5PARS=0 

SOLVE SET OF SIMULTANEOUS EQUATIONS 

226 CALL AMATF ( NCOL , I END , MK , KSP ARS , UMATRX ) 

227 CALL TOPLT ( T5TA , NT , 5STA , NQ , SFCN ) 

C WRITE (6,457) VALUA(LA) 

WRITE (9,465) ( S FCN ( L ) , L =1 , NQ ) 

WRITE (9,469) ( TFCN ( L ) , L =1 , NT ) 

Cl = 2.3025S5092 
C2 = 1. 0/(Cl>^ALPHL) 

C3 = 1 . 0/(Cl>^ALPHH) 

C4 □ C2^C2 
C5 — C3^C3 

UNKN(l) = SFCN(1)-C4^SFCN(3) 

UNKN(2) = SFCN(2)-C2^SFCN(3) 

UNKN(3) □ C4XSFCN(3) 

UNKN(4) = SFCN(1)-C5^SFCN(4) 

UNKN(5) = SFCN(2)-C3XSFCN(4) 

UNKN(6) = C5^SFCN(4) 

UNKN(7) = UNKN(2)+C1^ALPHL^UNKN(3) 

UNKN(8) □ UNKN(3)^C1^C1^ALPHL^ALPHL 
UNKN(9) = UNKN(5)+C1^ALPHH>^UNKN(6) 

UNKN(IO) = UNKN(6)^C1^C1)^ALPHHXALPHH 
C WRITE (7,524) UNKN 

5DEL = (SSTA(NS)-5STA(l)+0.2)/50.0 

C ADDING .1 TO GIVE SOME EXTRAPOLATION 18 NOV 76 
STRL(l) = 5STA(l)-*0.2 
DO 231 K = 2,50 
STRL(K) □ STRL(K-l) + SDEL 
231 CONTINUE 
KPLT = 0 
DO 232 J □ 1,50 
ALSTRL(J) = 10 . 0XXSTRL( J) 

SVAL = 10.0>^^(STRL(J)-SSTA(2)) 

SVLG = DLOGIO(SVAL) 

GLOW( J) = SFCN(1)+SFCN(2)KSVLG+SFCN(3)^(C4^SVAL^^ALPHL-C4-C2^SVLG) 
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GHIGHC J) = SFCNCl )+SFCN(2))^'SVLG+SFCN(4))^(C5X5VAL)<^ALPHH-C5-C3^SVLG) 
IF CSTRU J) .GT.S5TAC2) ) GO TO 229 

GPLTC J) = SFCN(1) + 5FCN(2)XSVLG+5FCN(3)^(C‘+^5VAL^^ALPHL-C^-C2^SVLG) 
GO TO 228 

229 GPLT(J) = SFCN(1)+SFCN(2)^SVLG+SFCN(4)>«(C5^SVAL^>^ALPHH'-C5-“C3^SVLG) 
228 KPLT = KPLT + 1 

232 CONTINUE 
KPLOT = 0 

DO 233 J = l,NODT 

5VAL = 10.0^^(WALS(J)-SSTA(2)) 

SVLG = DLOGIO(SVAL) 

IF (WALS(J) .GT.SSTAC2)) GO TO 23^ 

GPLATC J) = 5FCN(1)+SFCN(2)^SVLG + 5FCN(3)>((CAHSVAL^^ALPHL-CA-C2^5VLG) 
GO TO 235 

23A GPLAT(J) □ SFCN ( 1 ) +SFCN ( 2 ) S VL G + SFCN ( <^ ) « ( C5«S VAL ALPHH-C5-C3^S VLG ) 
235 KPLOT = KPLOT + 1 

233 CONTINUE 

PLOT THE ISOTHERMALS OF INTEREST THAT ARE NOT SAME AS STATIONS 

IF (NISO. EQ. NT) GO TO 21A 
DO 211 ISO=l,NISO 

PISO(ISO) □ D+E^CISOT(ISO)+F^CISOT(ISO)^^2 
KP = 0 

TRIPC ISO) = ( (CISOK ISO)-TSTA(NP) )/(T5TA(NP)+A60 . ) )^^2 
VALUMC ISO) = VALUA( LA))<( 1 . 0 - AKON^TRIP ( I SO ) ) 

DO 215 15=1 .50 

IF C IS . GT . Kf- LT ) GO TO 213 

ALGTCCIS) = (GPLT(IS)-P 150(1 SO ))/(!. 0 + VA L UM ( I SO ) x P I SO ( 1 5 0 ) ) 

IF (ALGTCCIS) .LT. 1.0. OR. ALGTC(IS).GT. 5.0) GO TO 213 
KP = KP+1 

ALSPE (KP) = STRL(IS) 

ALTFE (KP) = ALGTC (IS) 

206 k'PITE ( 7,507 ) V A LU A ( L A ) , Cl SOT ( ISO ) , ALTT ( IS ) , STRL ( I S ) , ANLO ( IS ) 

213 CONTINUE 

IVGRC2) □ KP 
IVGR(3) = 2 
IVGR(A) = 0 

CALL GPLOT ( ALTPE , AL SP E , I VGR ) 

211 CONTINUE 

PLOT THE PREDICTED ISOTHERMALS— LOG STRESS VS LOG TIME 

21A PCUT = .05 
IVGR(2) = 1 
IVGR(3) □ 3 
IVGR(A) □ 186 

CALL GPLOT ( CUT L OG , PCUT , I VGR ) 

DO 265 K=1,NT 
KP = 0 

TRIP(K)=( (TSTA(K)-TSTA(NP) )/(TSTA(NP)+A60 . ) )^^2 
VALUM(K)=VALUA(LA)^(1 . 0- AKON^TRIP ( K ) ) 

DO 260 1=1,50 
IF (I.GT.KPLT) GO TO 260 

ALGTC(I) = (GPLT(I)-‘TFCN(K))/(1.0 + VALUM(K)^TFCN(K)) 

IF (ALGTC(I).LT.1.0.0R.ALGTC(I).GT.5.0) GO TO 260 
KP=KP+1 

ALOGP(KP)=STRL(I) 

ALTP(KP)=ALGTC(I) 

C262 WRITE (7,507) VALUA ( L A ) , TSTA ( K ) , ALGTCC I ) , STRL ( I ) , TFCN ( K ) 

260 CONTINUE 

CALL INTENS(40) 

IVGR(2) = KP 
IVGR(3) = 0 

CALL GPLOT ( A L TP , ALOGP , I VGR ) 

CALL INTENS(2) 

265 CONTINUE 

C CALL ENDID(1,0,GRNAM) 

CALL DISPLA(l) 
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END OF ISOTHERMAL PLOT 


PLOT STRESS VS. LOG TIME FOR SELECTED TEMPERATURES 
SINCE THERE ARE NO PREDICTIONS, USE ONLY A=0 

256 IF (0PTC2)) GO TO 261 

IF (VALUA(LA).NE.O.O) GO TO 261 
CALL BEGIDC2) 

NUMPT = NODTW 
GRAXVC3) = 0. 

GRAXV(^) = 0.5 
GRAXVC5) = ^.5 
GRAXVC6) = .5 
GRAXVC7) = 2. 

GRAXVC8) - -1. 

GRAXVC9) = 0. 

RMIN =0.0 
RMAX □ 12.5 

IF (SSTA(3) .GT.1.00) RMAX=50. 

IF (S5TA(3) .GT.1.70) RMAX=75. 

IF (S5TA(3) .GT.1.88) RMAX=100. 

IF (SSTA(3) .GT.2.0I) RMAX=125. 

IF (55TA(3) .GT.2.10) RMAX=150 . 

IF (5STA(3) .GT.2.1S) RMAX=200, 

GRAYV(3) = 90. 

GRAYV(A) = RMIN 
GRAYVC5) = RMAX 
GRAYVC6) = 9. 

CALL XAXIS ( .9, .6,GRAXV) 

CALL YAXIS ( .9, .6,GRAYV) 

IVGR(3) = 3 

CALL TITLE (1 , 6 , 25 , T I T EL ) 

CALL CHARS ( 1 6 > DAT ES , 0 . , 1 . 7 , 9 . A , 12 ) 

CALL TITLE ( 3 , 1 2 , 2 0 , Q L ABX8 ) 

CALL CHARS ( 8 , GRL EG5 > 0 . , 8 . 0 , 9 . 5 , 1 5 ) 

IA = 1 

DO 26^ K=l,NISO 
IK = 0 

DO 263 I=l,NODTW 

IF (WTMP(I) .NE.CISOT(K)) GO TO 263 
IK=IK+1 

237 QXPL(IK)=WALT(I) 

QYPL(IK)=SIG(I) 

263 CONTINUE 
C 

C PLOT THE DATA FOR THE K TH TEMPERATURE 
C 

IVGRC2) = IK 
IVGR(^) = KSYM(K) 

CALL GPLOT ( QXP L , QYPL , I VGR ) 

26^ CONTINUE 

CALL CHARS ( 8 , GRL EG6 , 0 . , 8 . 5 , 9 . 5 > 15 ) 

CALL TITLE ( ^ , 16 , 20 , Q L ABXl) 

IVGRC2) = 1 
IVGR(^) = 186 

CALL GPLOT ( CUTLOG, PCUT, IVGR) 

C CALL ENDID(2,0,GRNAM) 

CALL DISPLA(l) 

C 

C END OF RAW DATA PLOT 
C 

C SET UP A TABLE OF LOG G = LOG H VALUES FOR THE CASE OF A=0 ONLY 
C 

C WRITE (9,A81) C XFI LS ( I I) , I I =1 , 49 , 3 ) 

C WRITE (9,481) ( YFI L S ( I I) , I I =1 , 50 , 3 ) 
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261 IF (VALUA(LA) ,NE. 0 . 0) GO TO 251 
DO 266 1=1,50 
DAZERO(I) = GPLAT(I) 

AMATRX(I,LA) = DAZERO(I) 

266 CONTINUE 
GO TO 252 

C 

C SET UP A TABLE OF LOG S AND LOG H VALUES FOR THE OTHER A VALUES 
C 

251 DO 267 1=1,50 

AMATRX (I, LA) = GPLAT(I) 

BNONZ(I) = GPLATCI)“ALTT(I) 

BNATRX (I, LA) = BNONZCI) 

267 CONTINUE 
C 

C PLOT THE G STATION FUNCTION VS LOG STRESS 
C 

252 5PLPT = SSTAC2) 

5PLY = 5FCNC1) 

SPLX = 10 . 0^>fS5TA(2) 

IF (0PTC3)) GO TO 254 
CALL BEGIDC3) 

NUMPT = KPLT 

CALL SCLBAK CXAXV , NUMPT , STRL , XMIS ) 

CALL GINTVL ( XNI 5 ( 1 ) , XMI S ( 2 ) , 1 0 , 0 , XRMIN , XRMAX ) 

GRAXVC3) = 0. 

GRAXVC4) = XRniN 
GRAXVC5) □ XRHAX 
GRAXVC6) = 10. 

CALL SCLBAK C YAXV , NUMPT , GP L T , YMI 5 ) 

CALL GINTVL ( YM I S ( 1 ) , YMI S ( 2 ) , 1 0 , 1 , YRMI N , YRMAX ) 

GRAYVC3) = 90. 

GRAYV(A) = YRMIN 
GRAYVC5) = YRMAX 
GRAYVC6) = 10. 

CALL XAXIS ( . 9, .6,GRAXV) 

CALL YAXI5 C . 9, .6,GRAYV) 

IVGRC2) □ KPLT 
IVGRC3) □ 0 
IVGRC4) = 70 

CALL GPLOT ( STR L , GP L T , I VGR ) 

C U'RITE (6,519) ( ( STR L ( I ) , I = 1 , 6 ) , ( GP L T ( I ) , I = 1 , 6 ) ) 

IVGRC2) □ 1 
IVGRC3) = 3 
IVGRC4) □ 71 

257 CALL GPLOT ( SP L PT , SP L Y , I VGR ) 

CALL TITLE ( 1 , 6 4 , 25 , T I T EL ) 

CALL CHARS ( 1 6 , DA T ES , 0 . , 1 . 7 , 9 . 4 , 1 2 ) 

CALL TITLE ( 3 , 4 , 20 , QL ABY2 ) 

CALL CHARS C 8 , GR L EG5 , 0 . , 8 . 0 , 9 . 5 , 1 5 ) 

CALL NUMBER ( 4 , V A L U A ( L A ) , 8 , 2 , GR L EG6 ) 

CALL CHARS ( 8 , GR L EG6 , 0 . , 8 . 5 , 9 . 5 , 1 5 ) 

CALL TITLE ( 4 , 1 6 , 2 0 , Q L AB Y1 ) 

C CALL ENDID(3,0,GRNAM) 

CALL DISPLA(l) 

C 

C PLOT THE G STATION FUNCTION VS STRESS 
C 

254 IF (OPT(4)) GO TO 253 
CALL BEGID(4) 

NUMPT □ KPLT 

CALL SCLBAK ( XAXV , NUMPT , A L STR L , XMI 5 ) 

CALL GINTVL ( XMI S ( 1 ) , XMI S ( 2 ) , 1 0 , 0 , XRMI N , XRMAX ) 

GRAXVC3) = 0. 

GRAXV(4) = XRMIN 
GRAXVC5) = XRMAX 
GRAXVC6) = 10. 

CALL SCLBAK ( YAXV , NUMPT , GPLT , YMIS ) 

CALL GINTVL ( YMI S ( 1 ), YMIS ( 2 ), 1 0 , 1 , YRMIN , YRMAX) 
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GRAYVC3) = 90. 

GRAYV(A) = YRMIN 
GRAYVC5) = YRMAX 
GRAYVC6) = 10. 

CALL XAXI5 ( .9, .6,GRAXV) 

CALL YAXI5 ( .9, .6,GRAYV) 

IVGRC2) = KPLT 
IVGRC3) = 0 
IVGR(A) = 70 

255 CALL GPLOT ( ALSTRL , GPLT , IVGR) 

C WRITE (6,519) ( ( AL5TRL ( I) , I = 1 , 6 ) , ( GPLT ( I ) , I = 1 , 6 ) ) 

IVGR(2) = 1 
IVGRC3) = 3 
IVGR(A) = 71 

CALL GPLOT ( 5PLX , 5PLY , I VGR ) 

CALL TITLE ( I , 6A , 25 , TITEL) 

CALL CHARS ( 16 , DAT E5 , 0 . , 1 . 7 , 9 . 4 , 12 ) 

CALL TITLE ( 3 , , 2 0 , Q L ABY2 ) 

CALL CHARS C 8 , GRL EG5 , 0 . , S . 0 , 9 . 5 , 1 5 ) 

CALL NUMBER ( A , V AL U A ( L A ) , 8 , 2 , GRL EG6 ) 

CALL CHARS ( 8 , GR L EG6 , 0 . , 8 . 5 , 9 . 5 , 1 5 ) 

CALL TITLE ( A , 12 , 20 , QL ABX8 ) 

C CALL ENDIDCA, 0,GRNAM) 

CALL DISPLA(l) 

PLOT THE P STATION FUNCTION VS TEMP 

253 TFIL(I) = TSTA(l) 

DTFIL(l) = TSTACI) 

PFIL(l) = TFCN(l) 

TRFIL(l) = 1000 ./(TSTA(1)+A60 . ) 

TDEL = ( TSTA(NT) - TSTACI) )/ 30.0 
249 NFIL = 31 

DO 269 1=2, NFIL 

TFIL(I) = TFIL(I-l) + TDEL 

DTFIL(I) = TFIL(I) 

TRFILCI) = 1000./(TFIL(I)+460.) 

248 PFIL(I) = D+E)^TFILCI) + F^TFILCI)^x2 
269 CONTINUE 

WRITE (7,487) ( T FI L ( I ) , I = 1 , N FI L ) 

WRITE (7,491) (PFIL(I), 1=1, NFIL) 

258 IF (0PT(5) ) GO TO 259 
CALL BEGIDC5) 

NUMPT = NFIL 

CALL SCLBAK (XAXV , NUMPT , TFI L , XMIS ) 

CALL GINTVL ( XMIS ( 1 ) , XMIS ( 2 ) , 1 0 , 0 , XRMIN , XRMAX ) 
GRAXVC3) = 0. 

GRAXVC4) = XRMIN 
GRAXVC5) = XRMAX 
GRAXVC6) = 10. 

CALL SCLBAK ( YAXV , NUMPT , PFI L , YMIS ) 

CALL GINTVL (YMISC 1 ) ,YMIS(2) , 10 , 1 , YRMIN, YRMAX) 
GRAYV(3) = 90. 

GRAYVC4) = YRMIN 
GRAYVC5) = YRMAX 
GRAYVC6) = 10. 

CALL XAXIS ( .9, .6,GRAXV) 

CALL YAXIS ( .9, .6,GRAYV) 

IVGRC2) = NFIL 
IVGRC3) = 0 
IVGRC4) = 70 

CALL GPLOT ( TFIL , PFIL , IVGR) 

IVGRC2) = NT 
IVGRC3) = 3 
IVGRC4) = 70 

CALL GPLOT ( QXPLTT , QYPLTT , IVGR) 

CALL TITLE ( 1 , 64 , 25 , T ITEL ) 

CALL CHARS ( 16 , DATES , 0 . , 1 . 7 , 9 . 4 , 12 ) 

CALL TITLE ( 3 , 4 , 20 , QL ABY3 ) 
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CALL CHARS ( 8 , GR L EG5 , 0 . , 8 . 0 , 9 , 5 , 1 5 ) 
CALL NUMBER ( ^ , V A L UA ( L A ) , 8 , 2 , GR L EG6 ) 
CALL CHARS ( 8 , GRL EG6 , 0 , , 8 . 5 , 9 . 5 . 15 ) 
CALL TITLE ( ^ ^ , 2 0 , QL ABX3 ) 

CALL ENDIDC5, O^GRNAM) 

CALL DISPLA(l) 


PLOT THE P STATION FUNCTION VS RECIPROCAL TEMP 


9 IF (0PTC6)) GO TO 268 
CALL BEGIDC6) 

NUMPT = NFIL 

CALL SCLBAK ( XAXV , NUMPT , TR FI L , XMI S ) 

CALL GINTVL C XMI S ( 1) , XMI S ( 2 ) , 1 0 , 0 , XRMIN , XRMAX ) 
GRAXV(3) = 0. 

GRAXV(A) = XRMIN 
GRAXVC5) □ XRMAX 
GRAXVC6) = 10. 

CALL SCLBAK ( YAXV , NUMPT , P FI L , YMI S ) 

CALL GINTVL ( YMI S ( 1 ) , YMI S ( 2 ) , 1 0 , 1 , YRMI N , YRMAX ) 
GRAYVC3) = 90. 

GRAYV(A) = YRMIN 
GRAYVC5) = YRMAX 
GRAYVC6) = 10. 

CALL XAXIS ( .9, .6,GRAXV) 

CALL YAXIS ( . 9, .6>GRAYV) 

IVGRC2) = NFIL 
IVGRC3) = 0 
IVGR(A) = 70 

CALL GPLOT ( TR FI L , P FI L , I VGR ) 

IVGR(2) = NT 
IVGRC3) = 3 
IVGR(A) □ 70 

CALL GPLOT ( QXP L RT , QYP L TT , I VGR ) 

CALL TITLE ( 1 , 6 A , 2 5 , T I T EL) 

CALL CHARS ( 16 , DATES , 0 . .. 1 . 7 , 9 . A , 12 ) 

CALL TITLE ( 3 , A , 2 0 , Q L ABY3 ) 

CALL CHARS ( 8 , GR L EG5 , 0 . , 8 . 0 , 9 . 5 , 1 5 ) 

CALL NUMBER ( A , V A L U A ( L A ) , 8 , 2 , GR L EG6 ) 

CALL CilARS (8, GRLEG6,0. ,8.5,9.5,15) 

CALL TITLE ( A , 28 , 20 , QLABXA ) 

CALL ENDIDC6 , 0 ,GRNAM) 

CALL DISPLACl) 

PLOT THE MASTER CURVE LOG SIGMA VS. THE PARAMETER G 


NOTA BENE MEG16B 

IF TEMPERATURE IS NOT A STATION VALUE, USE QUADRATIC CURVE FIT 
TO INTERPOLATE TEMP FUNCTION IN CALCULATION OF PARAMETER 


268 DO 285 I=l,NODT 
DO 270 IT=1,NT 
IY=IT 

IF (WTMP(I) .EQ.TSTA(IT)) GO TO 275 
270 CONTINUE 

PEST(I)=D+E)(UJTMP(I) + F^WTMP(I)K^2 
GO TO 280 

275 PEST( I ):=TFCN( lY) 

TRIPC I ) = ( (WTMPC I )-TSTA( HP) )/(TSTA(NP)+A6 0 . ) )>^^2 
VALUMC I ) =VALUA( LA)K( 1 . O-AKON^TRIP ( I ) ) 

28 0 GCALC(I)-IJALT(I) + PEST(I)>^(VALUM(I)>^U1ALT(I) + 1 . 0) 

285 CONTINUE 

286 IF (0PTC7)) GO TO 27A 
CALL BEGIDC7) 

NUMPT = KPLT 

CALL SCLBAK ( XAXV , NUMPT , GPLT , XMI5 ) 

CALL GINTVL ( XMI S ( 1 ) , XMI S ( 2 ) , 1 0 , 0 , XRMI N , XRMAX ) 
GRAXVC3) = 0. 
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GRAXV(^) = XRMIN 
GRAXVC5) = XRMAX 
GRAXVC6) = 10. 

CALL SCLBAK ( YAXV , NUMPT , STRL , YMI5 ) 

CALL GINTVL ( YMI5( 1 ) , YmS(2) , 10 , 1 , YRMIN, YRMAX) 

GRAYV(3) = 90. 

GRAYV(^) = YRMIN 
GRAYVC5) = YRMAX 
GRAYVC6) = 10. 

CALL XAXIS ( . 9, .6,GRAXV) 

CALL YAXIS ( . 9, .6,GRAYV) 

IVGRC2) □ 1 
IVGRC3) = 3 
IVGR(A) = 156 
SPLY = .6 

CALL GPLOT ( SPLY , SPLPT , I VGR ) 

IVGRC2) = KPLT 
IVGRC3) = 0 

CALL GPLOT ( GP LT , 5TRL , I VGR ) 

IVGRC2) = NODT 
IVGRC3) = 3 
IVGR(A) = 62 

CALL GPLOT ( GCALC , QALS , I VGR) 

CALL TITLE ( 1 , 6 A , 25 , TITEL ) 

CALL CHARS ( 1 6 , DATES , 0 . , 1 . 7 , 9 . A , 12 ) 

CALL TITLE ( 3 , 1 6 , 2 0 , QL ABYl ) 

CALL CHARS ( 8 , GR L EG5 , 0 . , 8 . 0 , 9 . 5 , 1 5 ) 

CALL NUMBER ( A , VALUA ( L A ) , 8 , 2 , GRL EG6 ) 

CALL CHARS ( 8 , GRL EG6 , 0 . , 8 . 5 , 9 . 5 , 15 ) 

CALL TITLE ( A , 28 , 20 , QL ABX5 ) 

CALL NUMBER ( 2 , S FCN ( A ) , 1 2 , 3 , GRL EGl ) 

CALL CHARS ( 1 2 , GR L EGl , 0 . , 1 . 1 , 2 . 0 , 1 5 ) 

CALL NUMBER ( 2 , S FCN ( 3 ) , 12 , 3 , GRL EGl) 

CALL CHARS ( 1 2 , GR L EGl , 0 . , 1 . 1 , 2 . A , 1 5 ) 

CALL NUMBER ( 2 , SFCN ( 2 ) , 12 , 3 , GRL EGl) 

CALL CHARS ( 1 2 , GRL EGl , 0 . , 1 . 1 . 2 . 8 , 15 ) 

CALL NUMBER C 2 , SFCN ( 1) , 12 , 3 , GRL EGl) 

CALL CHARS ( 1 2 , GR L EGl , 0 . , 1 . 1 . 3 . 2 , 1 5 ) 

CALL NUMBER ( 2 , QYPLTT( NT ) . 12, 3 . GRLEGl ) 

CALL CHARS ( 1 2 , GR L EGl , 0 . , 1 . 1 , 3 . 6 , 1 5 ) 

CALL NUMBER ( 2 , Q YP L TT ( 2 ) , 1 2 , 3 , GRL EGl ) 

CALL CHARS C 12 , GRL EGl , 0 . , 1 . 1 , A . 0 , 15 ) 

CALL NUMBER ( 2 , QYPLTTC 1 ) , 12 , 3 , GRLEGl ) 

CALL CHARS ( 1 2 , GRL EGl , 0 . , 1 . 1 , A . A , 15 ) 

C CALL ENDID(7,0,GRNAM) 

CALL DISPLA(l) 

C 

27A IF (0PT(8)) GO TO 281 
C 

281 IF (OPT(9)) GO TO 277 
C 

C CALCULATE STANDARD DEVIATION OF DATA USED— IN TERMS OF LOG TIME 
C 

277 CONTINUE 
SUMSQ=0 . 0 
SUMRES=0.0 
SUMAR=0 . 0 
SUMSQR =0.0 
SUMZD =0.0 
SUM2D2 =0.0 
JR = 1 

WRITE (9,A59) VALUA(LA) 

WRITE (9,A75) 

DO 325 J=l,NODUW 
YVAL = WALS(J) 

288 CALL TUSTR ( SSTA , YVAL , NQ , ALPHL , ALPHH, COEFGG) 

C WRITE (7,A50) ( COEFGGC KK) , KK=1 , NS ) 

YVALU = WTMP(J) 

321 CALL STACN ( TSTA , YVALU , NT , COEFP , N ) 


(STDEV) 
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C WRITE (7,<+5 0) (COEFP(KK) ,KK = 1,NT) 

PTRM=C0EFP(N-1 )^TFCN(N-l)+C0EFP(N))^TFCN(N)+C0EFP(N+l)KTFCNCN+l)+C0EFP- 
l(N+2)KTFCN(N+2) 

TRIPC J )=( (WTMPC J)-TSTA(NP) )/(T5TA(NP)+^60 . ))^^2 
VALUMC J)=VALUA(LA)X(1. O-AKON^TRIP(J)) 

H0UR5C J )=(GPLAT( J )”PTRM)/(1 , O+VALUMC J ) ^PTRM) 

RE5( J)=WALT( J )”H0UR5( J) 

SUNSQ = SUM5Q + RES( J)>^RE5( J) 

SUnRE5=5UMRE5+RES( J) 

ABRESC J )=DAB5(WALT( J)“H0URS( J) ) 

SUMAR-SUMAR+ABRES( J ) 

ZDEV(J) □ WALT(J)?<(1.0+VALUM(J)^PTRM)+PTRri-GPLAT(J) 

SUMZD = SUnZD + ZDEV(J) 

5UMZD2 = SUMZD2 + ZDEV ( J ) ^ZDEVC J ) 

WRITE (9,A89) WA L T ( J ), HOURS ( J RES ( J ), 5UMSQ , WTMP C J ), WAL5 ( J ), GPL AT ( J ), PTRM , N 
325 CONTINUE 

NOTA BENE MEG16B 

NO. OF PARAMETERS FITTED = ^+NT-1 SINCE TFCN AT ONE TEMP=0 

NO. OF DEGREES OF FREEDOM = NO. OF DATA - (4+NT-l) - 1 

NPAR = NT+NQ-1 
NDIV^^NODUW-NPAR-l 
IF (NDIV.LT.l) NDIV=1 
DOF=NDIV 

STDEV=SQRT(5UM5Q/D0F) 

SDZD=SQRT(SUMZD2/D0F) 

PSDZD(LA) □ SDZD^2.0 
TOT=NODUU 
AVDEV-SUMAR/TOT 
PSTDEVCLA)-STDEV 

526 WRITE (6,A85) S T DEV , DOF , SUMRES , A VDEV 

CALCULATE RESIDUALS AND NORMAL DEVIATES FOR ALL DATA (RSDEV) 

WRITE (9,A83) TITEL 
WRITE (9,A57) VALUA(LA) 

WRITE (9,A85) STDEV , DOF, SUMRES , AVDEV 

SMSQL^O , 0 

SMARSL^O . 0 

DO 375 J=l,NODTW 

YVALU = WTMP(J) 

351 CALL 5TACN ( T S T A , YV A L U , NT , CO EFP , N ) 

PTRM = COEFP(N-1)kTFCN(N- 1 )+C0EFP(N))^TFCN(N)+C0EFP(N + l))^TFCN(N + l)+C0EFP- 

l(M+2)>^TFCN(N+2) 

TRIPC J ) = ( (WTMPC J )-TSTA(NP) )/CTSTA(NP) + A60 . ) )^^2 
VALUMC J)=VALUA(LA)^(1 . 0 - AKON ^ TR I P ( J ) ) 

HOURC J )=(GPLAT( J )~PTRM)/(1 . O+VALUMC J)^PTRM) 

RESD( J )nWALT( J )-HOUR( J) 

IF (LA.NE.l) GO TO 35A 
PLRESD(J) = RESD(J) 

35 A R5DEVC J )=RESD( J )/STDEV 
370 CONTINUE 
375 CONTINUE 

WRITE (9,A95) WTMP ( J ). WA L S ( J ), WA L T ( J ), HOUR ( J ), RESD ( J ), RSDEV ( J ) 

WRITE (6,500) STDL,DENL 
WRITE (6,505) AVRESL , AVARSL 

PLOT OBSERVED VERSUS PREDICTED VALUES OF LOG TIME TO FAILURE 

379 IF (OPT(IO)) GO TO 385 
CALL BEGID(IO) 

GRAXV(A) □ 0.5 
GRAXV(5) = 5.5 
GRAXV(6) = 10. 

GRAYV(A) = .5 
GRAYV(5) = 5.5 
GRAYV(6) = 10. 
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380 


C 


CALL XAXI5 (.9, .6,GRAXV) 

CALL YAXIS ( .9, .6,GRAYV) 

IVGRC2) = NODUW 
1VGRC3) = 3 
IVGR(9) = 63 

CALL GPLOT ( HOURS , Q AL T , I VGR ) 

IVGRC2) = 2 
IVGRC3) = 0 

CALL GPLOT ( QPX . QPY , I VGR ) 

CALL GPLOT ( QPXL . QPYL , I VGR ) 

CALL GPLOT ( QPXR , QPYR , I VGR ) 

CALL TITLE ( 1 , 6 9 , 25 , TITEL) 

CALL CHARS ( 16 , DAT ES , 0 . , 1 . 7 , 9 . A , 12 ) 
CALL TITLE ( 3 , 9 , 20 , QL ABY7 ) 

CALL CHARS ( 8 , GRL EG5 , 0 . , 8 . 0 , 9 . 5 , 15 ) 
CALL NUMBER ( A , VAL U A ( L A ) , 8 , 2 , GRL EG6 ) 
CALL CHARS ( 8 , GRL EG6 , 0 . , 8 . 5 , 9 . 5 , 15 ) 
CALL TITLE ( A , , 2 0 , QL ABX7 ) 

CALL CHARS ( 12 , QSDKEY , 0 . , 7 . , 1 . A , 15 ) 
CALL NUMBER ( A , STDEV , 8 , 2 > GRL EGl) 

CALL CHARS ( 8 , GRL EGl , 0 . , 8 . A , 1 . A , 15 ) 
CALL ENDID(10,0,GRNAM) 

CALL DISPLA(l) 


385 CONTINUE 

END OF THE LOOP FOR VARIOUS A VALUES 

THE NEXT FEW PLOTS WILL EVALUATE THE VARIOUS A VALUES USED 


381 

IF 

(OPT(ll) ) 

GO 

TO 

382 

382 

IF 

(0PT(12)) 

GO 

TO 

383 

383 

IF 

(0PT(13) ) 

GO 

TO 

38A 


PLOT THE A SCAN -- S. D. VERSUS A VALUES CHOSEN 


38A 


396 


C 


C 


IF (OPT(IA)) GO TO 376 
DO 396 1=1, LAK 
QVALUA(I) = VALUKCI) 

CONTINUE 

WRITE (7,511) QVALUA 
WRITE (7,509) PSTDEV,PSTDL 
CALL BEGID(IA) 

GRAXV(A) = -.2 
GRAXV(5) □ .2 
GRAXV(6) □ 8. 

GRAYV(A) = 0.0 
GRAYV(5) = .8 
GRAYV(6) = 8. 

CALL XAXIS ( .9, .6,GRAXV) 

CALL YAXIS ( .9, .6,GRAYV) 

IVGR(2) = LAK 
IVGR(3) = 3 
IVGR(A) = 66 

CALL GPLOT ( QVALUA , PSTDEV , I VGR ) 

CALL TITLE ( 1 , 6 A , 25 , TITEL ) 

CALL CHARS ( 16 , DAT ES , 0 . , 1 . 7 , 9 . A , 12 ) 
CALL TITLE ( 3 , A , 2 0 , Q L ABYA ) 

CALL CHARS ( 8 , GRL EG5 , 0 . , 8 . 0 , 9 . 5 , 15 ) 
CALL NUMBER ( A , VALUA ( L A ) , 8 , 2 , GRL EG6 ) 
CALL CHARS ( 8 , GRL EG6 , 0 . , 8 . 5 , 9 . 5 , 15 ) 
CALL TITLE ( A , 8 , 20 , QL ABXA) 

CALL ENDID(1A,0,GRNAM) 

CALL DISPLA(l) 
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WRITE (7,523) ( PSDZDC J ) , J =1 , L AK ) 

WRITE (7,502) ( PSTDEV ( J ) , J = 1 , L AK ) 

PLOT THE RESIDUALS VERSUS TEMPERATURE 

376 IF (0PT(15)) GO TO 377 

CALL BEGID(15) 

NUMPT = NODTW 

CALL SCLBAK ( XAXV , NUMPT , QTMP , XMIS ) 

CALL GINTVL ( XMIS ( 1 ) , XMIS ( 2 ) , 1 0 , 0 , XRMIN , XRMAX) 
GRAXV(l) = 10. 

GRAXV(2) = -1. 

GRAXV(3) = 0. 

GRAXV(^) = XRMIN 
GRAXV(5) = XRMAX 
GRAXV(6) = .5 
GRAXV(7) = 2. 

GRAXV(8) = -1. 

GRAXV(9) = 0. 

CALL SCLBAK ( YAXV , NUMPT , P L RESD , YMI5 ) 

CALL GINTVL (YMIS( 1 ) , YMIS(2) , 10 , 1 ,YRMIN,YRMAX) 
GRAYV(l) = 10. 

GRAYVC2) = -1. 

GRAYVC3) = 90. 

GRAYV(^) = -1.0 
GRAYVC5) = 1.0 
GRAYVC6) = .5 
GRAYV(7) □ 2. 

GRAYV(S) □ -1 . 

GRAYV(9) = 0. 

CALL XAXIS ( .9, .6,GRAXV) 

CALL YAXIS ( .9, .6,GRAYV) 

IVGR(l) = 7 
IVGRC2) □ NUMPT 
IVGRC3) = 3 
IVGRCA) = 62 
IVGRC5) □ 1 
IVGRC6) = 15 
IVGRC7) = 0 

CALL GPLOT ( QTMP , PLRESD, I VGR ) 

CALL TITLE ( 1 , 6 ^ , 25 , T I T EL ) 

CALL CHARS ( 1 6 , DA T ES , 0 . , 1 . 7 , 9 . A , 1 2 ) 

CALL TITLE ( 3 , 20 , 20 . QL ABYA ) 

CALL CHARS C 8 , GR L EG5 , 0 . , 8 . 0 , 9 . 5 , 1 5 ) 

CALL NUMBER ( ^ , V A L U A C L A ) , 8 , 2 , GR L EG6 ) 

CALL CHARS ( 8 , GR L EG6 , 0 . , 8 . 5 , 9 . 5 , 1 5 ) 

CALL TITLE ( A , , 2 0 , Q L A BX3 ) 

CALL EMDID(15,0,GRNAM) 

CALL DISPLACl) 

PLOT THE RESIDUALS VERSUS LOG STRESS 

377 IF (0PT(16)) GO TO 378 
C 

CALL BEGID(16) 

NUMPT □ NODTW 

CALL SCLBAK ( XAXV , NUMPT , Q A L S , XMI 5 ) 

CALL GINTVL ( XMI S ( 1) , XMI S ( 2 ), 1 0 , 0 , XRMIN , XRMAX ) 
GRAXVC3) = 0. 

GRAXV(A) = XRMIN 
GRAXVC5) = XRMAX 
GRAXVC6) = .5 

CALL SCLBAK ( YAXV , NUMPT , P L R ESD , YMI S ) 

CALL GINTVL ( YMI S ( 1 ) , YMI S ( 2 ) , 1 0 , 1 , YRMIN , YRMAX) 
GRAYVC3) = 90. 

GRAYV(A) □ -1 . 0 
GRAYVC5) = 1.0 


Ij 
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GRAYVC6) □ .5 

CALL XAXIS ( .9, .6,GRAXV) 

CALL YAXIS f .9, .6,GRAYV) 

IVGRC2) •- NUMPT 
IVGRC3) = 3 
IVGR(A) = 62 

CALL GPLOT ( Q A L S , P L R ESD , I VGR ) 

CALL TITLE C 1 , 6 , 2 5 , T I T EL) 

CALL CHARS (1 6 , DAT E5 , 0 . , 1 - 7 , 9 . ^ , 12 ) 
CALL TITLE ( 3 , 20 , 20 , OL A3YA ) 

CALL CflARS (8, GRLEG5,0. ,8.0,9.5,15) 
CALL NUi:nCR (^, VALUA(LA) ,8,2,GRLEG6) 
CALL CHAPS (8, GRLEG6,0. ,8.5,9.5,15) 
CALL TITLE ( , 12 , 20 , QL ABX2 ) 

C CALL EHDIDC16,0,GRNAM) 

CALL DISPLA(l) 


C 

C PLOT THE RESIDUALS VERSUS LOG LIFE 
C 

378 IF (OPT(17)) GO TO 399 


C 


c 

c 

399 


C 

395 


CALL BEGID(17) 

NUMPT □ NODTW 

CALL SCLBAK ( XAXV , NUMPT , QALT , XMIS ) 

CALL GINTVL ( XMI S ( 1 ) , XMI S ( 2 ) , 1 0 , 0 , XRMI N , XRMAX ) 
GRAXV(3) = 0. 

GRAXVC9) = XRMIN 
GRAXVC5) □ XRMAX 
GRAXVC6) □ .5 

CALL SCLBAK ( YAXV , NUMPT , P L R ESD , YMI S ) 

CALL GINTVL ( YMI S ( 1 ) , YMI S ( 2 ) , 1 0 , 1 , YRMI N , YRMAX ) 
GRAYV(3) □ 90. 

GRAYVC9) n -1.0 
GRAYVC5) □ 1 . 0 
GRAYVC6) = .5 
CALL XAXIS ( . 9, .6,GRAXV) 

CALL YAXIS ( .9, .6,GRAYV) 

IVGRC2) = NUMPT 
IVGRC3) = 3 
IVGRC9) = 62 

CALL GPLOT ( Q AL T , PLRESD , I VGR ) 

CALL TITLE ( 1 , 6 9 , 25 , T I T EL ) 

CALL CHARS ( 1 6 , DA T E5 , 0 . , 1 . 7 , 9 . 9 , 1 2 ) 

CALL TITLE C 5 , 20 , 20 , QLABY9 ) 

CALL CHARS ( 8 , GR L EG5 , 0 . , 8 . 0 , 9 . 5 , 1 5 ) 

CALL NUMBER ( 9 , V A L U A ( L A ) , 8 , 2 , GRL EG6 ) 

CALL CHARS ( 8 , GR L EG6 , 0 . , 8 , 5 , 9 . 5 , 1 5 ) 

CALL TITLE ( 9 , 16 , 20 , QLABXl ) 

CALL ENDID(17,0,GRNAM) 

CALL DISPLA(l) 

IF (NXAV.NE.l) GO TO 395 
LAK □ LAK-NXAV 
CALL TERM 
GO TO 100 


C 

c 

c 

c 

C INPUT FORMATS 
C 


900 FORMAT (16F5.0) 


15 IF (K-6) 30,25,20 

20 ANS=Z^2 

AN5=ANSXX3 


GO TO 65 
25 ANS=2XZ 


ANS=2^ANS^X2 
GO TO 65 
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30 IF (K--^) ^5, ^0,35 

35 ANS=2^2 

AN5=AN5^ANS 
GO TO 65 
40 ANS=2?^2?f2 

GO TO 65 

45 IF (K-2) 60,55,50 

50 AN5=2H2 

GO TO 65 
55 ANS=2 

GO TO 65 
60 AM5=1. 

65 FUNC=AN5 

RETURN 
C 

70 FORMAT ( IHO , 5X , 72HTOO MANY COEFFICIENTS ASKED FOR - ONLY 8 ARE A“ 

IVAILABLE YOU CANNOT GET 13) 

END 

C 


402 FORMAT (1615) 

4 04 FORMAT ( 2 1 5 , F9 . 1 , F6 . 3 , 4X , 1 7 L 1 ) 
408 FORMAT ( F5 . 0 , F7 . 3 , F8 . 2 , F5 . 0 ) 
410 FORMAT (16A4) 

C 

C OUTPUT FORMATS 
C 


455 

457 

459 

461 

463 
465 
46 9 
46 7 
471 

475 

477 

481 

483 

485 


487 

489 

491 

495 

500 

501 

502 

503 

504 

505 

506 

507 
509 
511 
513 

515 

516 

517 


FORMAT ( IHl , 5CX, 18A4 ) 

FORMAT (1H0,50X,22HTHE VALUE OF ’A’ = ,F7.3) 

FORMAT (lfiO,5X,25H DATA USED IN CALCULATION , 22X , 22HTHE VALUE OF »A' - 
1 - , F7 . 3) 

FORMAT (1H0,5X,27H LONG TIME DATA NOT USED ,20X,22HTHE VALUE OF ’A’ 
1 ^ . F7 . 3) 

FORMAT (1H0,1X,10HG ST AT I ON S/ 1 OX , 1 0 F 1 0 . 3/ ) 

FORMAT (1U0,1X,24H G COEFFICIENTS / 1 OX , 1 0 FI 0 . 3/ ) 

FORMAT (li;0,lX,24H P STATION FUNCTION /I OX , 1 0 FI 0 . 3/ ) 

FORMAT ( IHO , IX, lOliP STAT I ONS/1 OX , 1 0 FI 0 . 3 ) 

FORMAT ( IHO , 5X, 3h'NO . , 5X, 1 IHl EMPERATURE,9X,6HSTRESS,13X,4HTIME,9X,1- 
11}1L0G(STRFSS),9X,9I!L0G(TIME)/) 

FORMAT (lflO,20X,8HLOG TIME,22X,6H T EMP , 2X , 1 OHL OG STRESS , 

1 3X,5H G ,7X,3H P /11X,4H 0BS,6X,5H PRED,5X,5f! DI FF , 5X , 6HSUM SQ ) 
FORMAT (1HO,50HTHE HF,NP,NG,MK, N L , N T , N S , N I SO , NX A V VALUES ARE ,41- 
15,5X,515) 


F6 . 1,2X, 14H- 


,F7 . 3, 3X,26H- 


FORMAT ( IX, 18F7 . 3) 

FORMAT ( Ifll , COX, 18A4 ) 

FORMAT (1H0,16HS. D. OF RES - , F6 . 3 , 2X , 6 HDO F = 

ISl'M OF RES = , F6 . 3 , 2X, 14HAV ABS DEV □ ,F6.3) 

FORMAT riH0,7HrCMP = / 1 0 X , 1 0 F 1 0 . 3 ) 

FORMAT (8X,4F10.5,F10.1,3{-10.5,I4) 

FORMAT (1H0,71!P(T) “ / 1 0 X , 1 0 F 1 0 . 3 ) 

FORMA T (6X,F10.0,3F12.5,F12.4,ri2,3) 

FORMAT ( 1110 , 39I1STANDARD ERROR OF PREDICTED VALUES = 

INO. OF DATA FREDICTED - , F6 . 1 ) 

FORMAT ( IHO , 1 IFl 0 . 4 ) 

FORMAT (IX, *S. D. OF DATA USED = ’,9F8.3) 

^O^M^T (IX, *S. D. OF DATA PRED. = ’,9F8.3) 

FORMAT ( IHO , IX, ’CUTOFF = ’ , 2 FI 1 . 2 , 3X , ’ REDUCED DATA BEGIN AT LIFE = 
1 r]1.2,3X,F6.2,2X,F5,0) 

FC2.MAT (1X,39HAVG, DEVIATION OF PREDICTED VALUES ° , F7 . 3 , 3X , 39H- 

lAVG. ABS. DEV. OF PREDICTED VALUES = , F7 . 3 ) 

rOFMIAT (IHO, IX, ’NO. OF DATA USED = ’ , 1 6 , 5X , ’ NO . WEIGHTED =’, 

1 I6,5X,’NO. PREDICTED =’,I6,5X,’D 0 F =’,F6.0 > 

FORMAT (4X,F5.2, F10..0,3X,FS.3,3X,F10.4,2X,F10.4) 
rORMAT (5X, 9F6 . 3,2X, 9F5 . 3) 

FORMAT (2X,23HTIIE VALUES OF A = ,9F8.3) 

FORMAT (3X,54HA VALUE TEMP LOG TIME 

FORMAT ( IHl ,2X, 3CHRESULTS FROM MEG16B (MAY 78) 

1 15X.39HL0GH + A LOGH ^ P(T) + P(T) = G(SIGMA) 

FORMAT (10X,4A4) 

FORMAT (2X,12F6.2) 


LOG S 
/ 

) 


STRESS 


) 


29 


oooooooo oo ooo ooo ooo 


518 F0I:MAT (2X, 1216 ) 

519 FORMAT (2X, 12F8. 3) 

523 FORMAT (IX, D. OF ZERO DEV. = »,9F8.3) 

529 FORMAT (IX, 11F9.3) 

C 

399 STOP 
C 

END 

C 

c 

SUBROUTINE TUSTR (XS , X , NS , ALPL , ALPH , COEF) 

C 

C THIS PROGRAM SETS UP THE STRESS FUNCTION IN TERMS OF 
C TWO STRAIGHT LINES MERGED TOGETHER SMOOTHLY 
C 

C XS IS THE VECTOR OF ’STATION VALUES’ ONLY XS(2) IS USED IN THIS VERSION 

C X IS THE VALUE OF LOG STRESS TO BE EXPRESSED IN TERMS OF THE MODEL EQ. 

C NS IS THE NUMBER OF STATION VALUES SET = 3 USUALLY 
C ALPL IS THE EXPONENT FOR THE LOW STRESS REGION 

C ALPH IS THE EXPONENT FOR THE HIGH STRESS REGION 

C COEF IS THE VECTOR OF COEFFICIENTS FOR EACH X VALUE 
C 

DOUBLE PRECISION XS , X, COEF, S , SLOG, SSQ , Cl , C2 , C3 
DIMENSION XS(9) ,C0EF(9) 

C WRITE (7,200) XS 

10 S = 10 . 0^^KX-XS(2) ) 

SLOG □ DLOG10(10.0^>^X/10.0^XXS(2)) 

SSQ = S5<><2 
Cl a 2,302585093 
C2 = 1.0/(C15^ALPL) 

C3 = 1 . 0/(ClxALPH) 

IF (X.GT,XS(2) ) GO TO 20 

FOR STRESS LT MID STRESS REGION 1 

15 COEF( 1 ) □ 1 . 0 
C0EF(2) = SLOG 

C0EF(3) = C2^(C2>^5^^ALPL-C2-SL0G) 

COEFC9) = 0.0 
GO TO 30 

FOR STRESS GT MID STRESS REGION 2 

20 COEF(l) =1.0 
C0EF(2) = SLOG 
C0EF(3) = 0.0 

C0EF(9) = C3'><(C3>(S)<^ALPH-C3-SLOG) 

30 CONTINUE 

WRITE (7,209) X 
WRITE (7,201) COEF 
WRITE (7,203) S, SLOG, SSQ 

200 FORMAT (1X,3F8.9) 

201 FORMAT (10X,5F10.5) 

203 FORMAT (1X,9F9.5) 

209 FORMAT (lOX,* X = ’,F7.9) 

RETURN 
END 


SUBROUTINE STACN (X, Y, NS , COEF, N ) 

TO EXPRESS A SINGLE VALUE OF A VARIABLE IN TERMS OF DISCRETE ’STATIONS’ 
X IS THE VECTOR OF ’STATION VALUES’ 

Y IS THE VALUE TO BE EXPRESSED IN TERMS OF STATION VALUES 

NS IS THE NUMBER OF STATION VALUES 

COEF IS THE VECTOR OF COEFFICIENTS FOR EACH Y VALUE 
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C GIVEN A TABLE OF DISCRETE X VALUES AND ONE VALUE, CALLED Y , 

C FIND, FOR THE Y VALUE, THE COEFFICIENTS OF A TYPE OF LAGRANGIAN 
C INTERPOLATION EQUATION 
C 

DOUBLE PRECISION X , Y , CO EF , CA , CB , CC , CD , CE , CF , CH , C J , CK 
DIMENSION C0EF(32),X(32) 

WRITE (7,205) X 
WRITE (7,205) Y 
5 CA = 0.0 
CB □ 0.0 
CC □ 0.0 
CD = 0 . 0 
CE = 0.0 
CF = 0 . 0 
CH □ 0 . 0 
CJ = 0 . 0 
CK □ 0.0 
DO 10 K=l,32 
10 COEF(K) = 0.0 
30 DO 20 K= 1,NS 
M^NS-K 

WRITE (7,206) K,M 
IF (Y. LT.XC2) ) GO TO 15 
IF C Y. GE.XCNS-l ) ) GO TO 16 
IF (Y.GE.X(M) ) GO TO 17 
20 CONTINUE 

17 N = M 
CA = Y-X(N-l) 

CB □ Y-X(N) 

CC ° Y~X(N + 1 ) 

CD - Y-X(N+2) 

CF = X(N )-X(N+l) 

CE = X(N )-X(N-l) 

CH = X(N+1)-X(N-1) 

CJ - X(N )-X(N+2) 

CK - X(N+l)-X(N+2) 

^0 COEF(N-l) - ( CD^CC)/(CE^CH))^0 . 5 
COCFCN + 2) = (CB^^CC)/(CJ^CK)>^0.5 

COEFCN ) □ ( (CAHCC>^CJ) + (CC^CDHCE))/(CE><CF^CJ)*^0 .5 
C0I:F(N + 1) =~((CA^CBKCK) + (CB^CD><CH))/(CF^CH^CK)^0.5 
GO TO A5 

15 N-2 
GO TO 18 

16 N=MS-1 

18 COEF(N-l) = (X(N)-Y)>^CX(N + 1)-Y)/((X(N)-X(N-1))>^(X(N + 1)-X(N-1))) 

COEFCN) - (XCN + I )-Y)^ ( Y-XCN-l ) )/( CXCN + 1 )-X(N) )>< (X(N)-X(N-l ) ) ) 
COEFC M + 1 ) = (X( N)-Y)^( Y-X( N-1 ))/((X(N)-X(N + l))^(X(N + l)-X(N-l))) 

A5 CONTINUE 

WRITE (7,206) N 
WRITE (6,205) Y 

WRITE (6,205) ( CO EF ( L ) , L = 1 , N 5 ) 

205 FOR[*AT (IX, lOFlO.^) 

206 FORMAT (IX, 3IA) 

RETURN 

END 
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SUBROUTINE AMATF (M, N , MK , KSPARS , U ) 


C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 


5 

c 

c 

c 


10 

15 

C 

20 

C 

c 

c 

c 

c 

c 


25 

c 

c 

c 

c 

c 


30 

C 

C 

c 


MLSS/MFSS PROGRAM BEING USED TO SOLVE THE REDUNDANT EQS 

THE U MATRIX MUST BE BROKEN UP INTO THE PROPER MATRICES 

M IS THE NUMBER OF COLUMNS IN MATRIX 

N IS THE NUMBER OF ROWS 

MK IS THE COLUMN TO BE USED AS DEPENDENT VARIABLE 

KS IS A COUNTER NOT USED IN THIS VERSION 

U IS THE MATRIX TO BE INVERTED 

XF IS THE SOLUTION VECTOR 

DOUBLE PRECISION XF 

DOUBLE PRECISION F( 20 0 ) , AUT ( 2304 ) , TRAC( 20 0 ) , U ( 20 0 , 48 ) , A ( ^8 , 20 0 ) , XXT ( 48 , 48 ) 

COMMON /CCOMN/ XF(48),NP 

DO 5 1=1,48 

DO 5 J=l,200 

A(I,J)=0. 

READ X 

MKM1=MK-1 
MKP1=MK+1 
DO 20 1=1, N 
DO 10 J=1,MKM1 
A(J,I)=U(I, J) 

DO 15 J=MKP1,M 
AC J-1,I)=U(I, J) 

WRITE (6,60) (U(I, J), J=1,M) 

F(I)=“U(I,MK) 

CONTINUE 

M=M-1 


WRITE (6,60)((A(J,I),J=1,M),I=1,N) 
T 

FORM XX 

DO 25 1=1, M 
DO 25 J=1,M 
XXTCI, J)=0. 

DO 25 K=1,N 

XXTCI, J)=XXT(I, J)+A(I,K)XA(J,K) 

DO 7 1=1, N 

7 WRITE (6,60) (XXT( I , J ) , J=1 ,M) 

FORM XF 

DO 30 1=1, M 
XF(I)=0. 

DO 30 K=1,N 

XF(I)=XF(I) + A(I,K)>«F(K) 

PREPARE FOR LEAST SQUARES 
EPS=l.E-5 


K = 0 

DO 35 J=1,M 
DO 35 1=1, J 
K = K + 1 

35 AUT(K)=XXT(I, J) 

KTOT=K 

C WRITE (6,45) ( AUT ( K ) , K=1 , KTOT ) 

36 CALL DMFSS ( AUT , M , EPS , IRANK , TRAC) 
C WRITE (6,50) IRANK 

C WRITE (6,45) ( AUT ( K) , K=1 , KTOT) 
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COMPUTE 

T - 


A = 

(XX ) 

CALL 

DML55 ( 

:aut,m. 

WRITE 

(6,50) 

lER 

WRITE 

(6,45) 

(XF(I) 


1 

XF, THE SOLUTION OF 
IRANK,TRAC,0,XF,IER) 


XF 


T 

XX A= 0 


RETURN 

0 FORMAT (1H1/5H M= , I3,3X,AHN= , I3/1X/2H X/IX) 

5 FORMAT (2H F/1X/C8E16 .8) ) 

0 FORMAT CIS) 

0 FORMAT (2X,10F8.3) 

END 


SUBROUTINE DMFSS 
PURPOSE 

GIVEN A SYMMETRIC POSITIVE SEMI DEFINITE MATRIX , DMFSS WILL 

(1) DETERMINE THE RANK AND LINEARLY INDEPENDENT ROWS AND 
COLUMNS 

(2) FACTOR A SYMMETRIC SUBMATRIX OF MAXIMAL RANK 

(3) EXPRESS NONBASIC ROWS IN TERMS OF BASIC ONES, 

EXPRESS NONBASIC COLUMNS IN TERMS OF BASIC ONES 
EXPRESS BASIC VARIABLES IN TERMS OF FREE ONES 

SUBROUTINE DMFSS MAY BE USED AS A PREPARATORY STEP FOR THE 
CALCULATION OF THE LEAST SQUARES SOLUTION OF MINIMAL 
LENGTH OF A SYSTEM OF LINEAR EQUATIONS WITH SYMMETRIC 
POSITIVE SEMI-DEFINITE COEFFICIENT MATRIX 

USAGE 

CALL DMFSSCA.N, EPS, IRANK, TRAC) 

DESCRIPTION OF PARAMETERS 

A - UPPER TRIANGULAR PART OF GIVEN SYMMETRIC SEMI- 

DEFINITE MATRIX STORED COLUMNWISE IN COMPRESSED FORM 
ON RETURN A CONTAINS THE MATRIX T AND, IF IRANK IS 
LESS THAN N, THE MATRICES U AND TU 
A MUST BE OF DOUBLE PRECISION 
N - DIMENSION OF GIVEN MATRIX A 

EPS - TESTVALUE FOR ZERO AFFECTED BY ROUND-OFF NOISE 

IRANK - RESULTANT VARIABLE, CONTAINING THE RANK OF GIVEN 

MATRIX A IF A IS SEMI -0EFI N I T E 

IRANK = 0 MEANS A HAS NO POSITIVE DIAGONAL ELEMENT 

AND/OR EPS IS NOT ABSOLUTELY LESS THAN ONE 
IRANK =-l MEANS DIMENSION N IS NOT POSITIVE 
IRANK =-2 MEANS COMPLETE FAILURE, POSSIBLY DUE TO 
INADEQUATE RELATIVE TOLERANCE EPS 
TRAC - VECTOR OF DIMENSION N CONTAINING THE 

SOURCE INDEX OF THE I-TH PIVOT ROW IN ITS I-TH 
LOCATION, THIS MEANS THAT TRAC CONTAINS THE 
PRODUCT REPRESENTATION OF THE PERMUTATION WHICH 
IS APPLIED TO ROWS AND COLUMNS OF A IN TERMS OF 
TRANSPOSITIONS 

TRAC MUST BE OF DOUBLE PRECISION 

REMARKS 

EPS MUST BE ABSOLUTELY LESS THAN ONE. A SENSIBLE VALUE IS 
SOMEWHERE IN BETWEEN AND 10XJ<(-6) 

THE ABSOLUTE VALUE OF INPUT PARAMETER EPS IS USED AS 
RELATIVE TOLERANCE. 

IN ORDER TO PRESERVE SYMMETRY ONLY PIVOTING ALONG THE 
DIAGONAL IS BUILT IN. 

ALL PIVOTELEMENTS MUST BE GREATER THAN THE ABSOLUTE VALUE 
OF EPS TIMES ORIGINAL DIAGONAL ELEMENT 
OTHERWISE THEY ARE TREATED AS IF THEY WERE ZERO 
MATRIX A REMAINS UNCHANGED IF THE RESULTANT VALUE IRANK 
EQUALS ZERO 
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C 

C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 

C NONE 

C 

C METHOD 

C THE SQUARE ROOT METHOD WITH DIAGONAL PIVOTING IS USED FOR 

C CALCULATION OF THE RIGHT HAND TRIANGULAR FACTOR. 

C IN CASE OF AN ONLY SEMI-DEFINITE MATRIX THE SUBROUTINE 

C RETURNS THE IRANK X IRANK UPPER TRIANGULAR FACTOR T OF A 

C SUBMATRIX OF MAXIMAL RANK, THE IRANK X (N-IRANK) MATRIX U 

C AND THE (N-IRANK) X CN-IRANK) UPPER TRIANGULAR TU SUCH 

C THAT TRANSP05E(TU)xTU = I + TRANSP0SE(U)>^U 

C 

c 

c 

SUBROUTINE DMFSS (A , N , EPS , IRANK, TRAC) 

C 

C 

C DIMENSIONED DUMMY VARIABLES 

DIMENSION A(1),TRAC(1) 

DOUBLE PRECISION SUM, A , TRAC, PIV, HOLD 
C 

C TEST OF SPECIFIED DIMENSION 

IF(N)36.36,1 
C 

C INITIALIZE TRIANGULAR FACTORIZATION 

1 IRANK=0 
ISUB=0 
KPIV=0 

J = 0 

PIV=0.D0 

C 

C SEARCH FIRST PIVOT ELEMENT 

DO 3 K=1,N 

J=J + K 

TRAC(K)=A( J) 

IF(A(J)-PIV)3,3,2 

2 PIV=A(J) 

K5UB=J 

KPIV=K 

3 CONTINUE 
C 

C START LOOP OVER ALL ROWS OF A 

DO 32 1=1, N 

ISUB=ISUB+I 
IM1=I-1 
^ KMI=KPIV-I 
IF(KMI)35,9,5 
C 

C PERFORM PARTIAL COLUMN INTERCHANGE 

5 JI=KSUB-KMI 
IDC=JI-ISUB 
JJ=ISUB-IM1 

DO 6 K=JJ,ISUB 

KK=K+IDC 

HOLD=A(K) 

A(K)=A(KK) 

6 A(KK)=HOLD 
C 

C PERFORM PARTIAL ROW INTERCHANGE 

KK=KSUB 
DO 7 K=KPIV,N 
II=KK-KMI 
HOLD=A(KK) 

A(KK)=A(II) 

A(II)=HOLD 

7 KK=KK+K 
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PERFORM REMAINING INTERCHANGE 
JJ=KPIV-1 
II = ISUB 
DO 8 K=I,JJ 
H0LD=A(II) 

A(II)=A(JI) 

A( JI)=HOLD 
II=II+K 

8 JI=JI+1 

9 IF(IRANK)22,10,10 

RECORD INTERCHANGE IN TRANSPOSITION VECTOR 

10 TRAC(KPIV)=TRAC(I) 

TRAC(I)=KPIV 

MODIFY CURRENT PIVOT ROW 
KK=IM1-IRANK 
KMI=ISUB-KK 
PIV=0.D0 
IDC=IRANK+1 
JI=ISUB-1 
JK=KMI 
JJ=ISUB-I 
DO 19 K=I,N 
SUM=0.D0 

BUILD UP SCALAR PRODUCT IF NECESSARY 
IF(KK)13,13,11 

11 DO 12 J=KMI,JI 
SUM=SUM-A( J)^A( JK) 

12 JK=JK+1 

13 JJ=JJ+K 
IF(K-I)1^,1‘^>16 

1^ SUM=A(I5UB)+SUM 

TEST RADICAND FOR LOSS OF SIGNIFICANCE 
IF(SUM-DABS(A(ISUB)^DBLE(EP5) ) )20,20,15 

15 A(ISUB)=DSQRT(SUM) 

KPIV=I+1 

GOTO 19 

16 SUM=(A( JK)+SUM)/A(ISUB) 

AC JK)=SUM 

SEARCH FOR NEXT PIVOT ROW 
IF(A(JJ))19,19,17 

17 TRAC(K)=TRAC(K)-SUM^SUM 
HOLD=TRACCK)/A( JJ) 

IF(PIV-H0LD)18,19,19 

18 PIV=HOLD 
KPIV=K 
KSUB=JJ 

19 JK=JJ+IDC 
GOTO 32 

CALCULATE MATRIX OF DEPENDENCIES U 

20 IF(IRANK)21,21,37 

21 IRANK=-1 
GOTO ^ 

22 IRANK=IM1 
II=ISUB-IRANK 
JI=II 

DO 26 K=1,IRANK 

JI=JI-1 

JK=ISUB-1 

JJ=K-1 

DO 26 J=I,N 

IDC=IRANK 

SUM=0 .DO 
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Kru=Ji 

T. \ y 

IF( J j)25,25.23 
23 DO 2^ L=1 , JJ 
IDC=IDC-1 

SUri = 5UM-A(KMI))^A(KK) 
KMI=KMI-IDC 
2A KK:=KK-1 

25 A(KK)=(SUM+A(KK) )/A(KMI) 

26 JK=JK+J 
C 

C CALCULATE I + TRANSP05E(U)^^U 

JJ=I5UB-I 
PIV=0 . DO 
KK=I5UB-1 
DO 31 K=I,N 
JJ=JJ+K 
IDC^^O 

DO 28 J=K,N 
SUM-0 . DO 
KMI=JJ+IDC 
DO 27 L=II,KK 
JK=L+IDC 

27 SUM=SUM+A( L)^A( JK) 

A(Km )=SUM 

28 IDC=IDC+J 
A(JJ)=A( JJ)+1.D0 
TRAC(K)=A(JJ) 



SEARCH NEXT DIAGONAL ELEMENT 
IFCPIV-A(JJ))29,30,30 

29 KPIV'K 
KSUB=JJ 
PIV=A( JJ) 

30 II=II+K 
KK=KK+K 

31 CONTINUE 
GOTO A 

32 CONTINUE 

33 IF(IRANK)35,3A,35 
5A IRANK=N 

35 RETURN 


ERROR RETURNS 


RETURN IN CASE OF ILLEGAL DIMENSION 
36 IRANK=“1 
RETURN 


INSTABLE FACTORIZATION OF I+TRANSPOSE( U ) 

37 IRANK=-2 
RETURN 
END 

C 

c 

C SUBROUTINE DMLSS 

C 

C PURPOSE 

C SUBROUTINE DMLSS IS THE SECOND STEP IN THE PROCEDURE FOR 

C CALCULATING THE LEAST SQUARES SOLUTION OF MINIMAL LENGTH 

C OF A SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS WITH SYMMETRIC 

C POSITIVE SEMI-DEFINITE COEFFICIENT MATRIX, 

C 

C USAGE 

C CALL DMLSS(A,N, IRANK, TRAC, INC, RHS, lER) 

C 
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DESCRIPTION OF PARAMETERS 

A - COEFFICIENT MATRIX IN FACTORED FORM AS GENERATED 

BY SUBROUTINE MFSS FROM INITIALLY GIVEN SYMMETRIC 
COEFFICIENT MATRIX A STORED IN N^(N+l)/2 LOCATIONS 
A REMAINS UNCHANGED 
A MUST BE OF DOUBLE PRECISION 
N - DIMENSION OF COEFFICIENT MATRIX 

IRANK - RANK OF COEFFICIENT MATRIX, CALCULATED BY MEANS OF 
SUBROUTINE DMFSS 

TRAC - VECTOR OF DIMENSION N CONTAINING THE 

SUBSCRIPTS OF PIVOT ROWS AND COLUMNS, I.E. THE 
PRODUCT REPRESENTATION IN TRANSPOSITIONS OF THE 
PERMUTATION WHICH WAS APPLIED TO ROWS AND COLUMNS 
OF A IN THE FACTORIZATION PROCESS 
TRAC IS A RESULTANT ARRAY OF SUBROUTINE MFSS 
TRAC MUST BE OF DOUBLE PRECISION 
INC - INPUT VARIABLE WHICH SHOULD CONTAIN THE VALUE ZERO 
IF THE SYSTEM OF SIMULTANEOUS EQUATIONS IS KNOWN 
TO BE COMPATIBLE AND A NONZERO VALUE OTHERWISE 
RHS - VECTOR OF DIMENSION N CONTAINING THE RIGHT HAND SIDE 
ON RETURN RHS CONTAINS THE MINIMAL LENGTH SOLUTION 
RHS MUST BE OF DOUBLE PRECISION 
lER - RESULTANT ERROR PARAMETER 
lER □ 0 MEANS NO ERRORS 

lER =~1 MEANS N AND/OR IRANK IS NOT POSITIVE AND/OR 
IRANK IS GREATER THAN N 

lER = 1 MEANS THE FACTORIZATION CONTAINED IN A HAS 
ZERO DIVISORS AND/OR TRAC CONTAINS 
VALUES OUTSIDE THE FEASIBLE RANGE 1 UP TO N 

REMARKS 

THE MINIMAL LENGTH SOLUTION IS PRODUCED IN THE STORAGE 
LOCATIONS OCCUPIED BY THE RIGHT HAND SIDE. 

SUBROUTINE DMLSS DOES TAKE CARE OF THE PERMUTATION 
WHICH WAS APPLIED TO ROWS AND COLUMNS OF A. 

OPERATION IS BYPASSED IN CASE OF A NON POSITIVE VALUE 
OF IRANK 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NONE 

METHOD 

LET T, U, TU BE THE COMPONENTS OF THE FACTORIZATION OF A, 

AND LET THE RIGHT HAND SIDE BE PARTITIONED INTO A FIRST 
PART XI OF DIMENSION IRANK AND A SECOr^D PART X2 OF DIiMENSION 
N-IRANK. THEN THE FOLLOWING OPERATIONS ARE APPLIED IN 
SEQUENCE 

(1) INTERCHANGE RIGHT HAND SIDE 

(2) XI □ XI + U ^ X2 

(3) X2 =-TRANSrOSE(U) x XI 

(A) X2 □ INVERSE(TU) x INVERSEC TRANSPOSEC TU ) ) ^ X2 

(5) XI □ XI + U ^ X2 

(6) XI = INVERSE(T) K INVERSEC TRANSPOSEC T ) ) ^ XI 

(7) X2 =-TRAN5P0SE(U) XI 

(8) X2 = INVERSE(TU) ^ INVERSEC TRANSPOSEC TU ) ) ^ X2 
C9) XI = XI + U ^ X2 

(10)X2 n TRAN5P05ECU) X XI 
cm REINTERCHANGE CALCULATED SOLUTION 

IF THE SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS IS SPECIFIED 
TO BE COMPATIBLE THEN STEPS C2), C3), CA) AND C5) ARE 
CANCELLED. 

IF THE COEFFICIENT MATRIX HAS RANK N, THEN THE ONLY STEPS 
PERFORMED ARE Cl), C6) AND Cll). 
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SUBROUTINE DML SS ( A , N , IRANK , TRAC , INC , RHS , I ER ) 

C 

C 

C DIMENSIONED DUMMY VARIABLES 

DIMENSION A(1),TRAC(1),RHS(1) 

DOUBLE PRECISION SUM , A , RHS , TRAC , HOLD 
C 

C TEST OF SPECIFIED DIMENSIONS 

IDEF=N-IRANK 
IF(N)33,33,1 

1 IF(IRANK)33,33,2 

2 IF( IDEF)33,3,3 
C 

C CALCULATE AUXILIARY VALUES 

3 ITE = IRANK^( IRANK + D/2 
IX2=IRANK+1 

NP1=N+1 
IER = 0 
C 

C INTERCHANGE RIGHT HAND SIDE 

JJ = 1 
11 = 1 

^ DO 6 1=1, N 
J=TRAC(II) 

IF( J )31 , 31 , 5 

5 H0LD=RH5(II) 

RHS(II)=RHS(J) 

RHS( J)=HOLD 

6 II=1I+JJ 
IF(JJ)32,7,7 

C 

C PERFORM STEP 2 IF NECESSARY 

7 ISW=1 

IF(INC>(IDEF)8,28,8 

C 

C CALCULATE XI = XI + U ^ X2 

8 ISTA=ITE 

DO 10 I=1,IRANK 

ISTA=ISTA+1 

JJ=ISTA 

SUM=0 .DO 

DO 9 J=IX2,N 

SUM = SUM + A( JJ))^RHS( J) 

9 JJ=JJ+J 

10 RHS(I)=RHS(I)+SUM 
G0T0(11,28,11),1SW 

C 

C CALCULATE X2 = TRANSPOSECU) ^ XI 

11 ISTA=ITE 

DO 15 1=1X2, N 

JJ=ISTA 

SUM=0 .DO 

DO 12 J=1,IRANK 

JJ=JJ+1 

12 SUM=SUM+A(JJ)^RHSCJ) 

G0T0(13,13,1<^),ISW 

13 SUM=“SUM 

14 RHS(I)=SUM 

15 ISTA=ISTA+I 
GOTO(16,29,30),ISW 

C 

C INITIALIZE STEP (4) OR STEP (8) 

16 ISTA=IX2 
IEND=N 
JJ=ITE+ISTA 

C 

C DIVISION OF XI BY TRANSPOSE OF TRIANGULAR MATRIX 

17 SUM=0.D0 

DO 20 I=ISTA,IEND 
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IF(A(JJ))18,31,18 

18 RHS(I)=(RHS(I)-SUM)/ACJJ) 

IF(I-IEND)19,21,21 

19 JJ=JJ+ISTA 
SUM-0, DO 

DO 20 J=ISTArI 
5UM=5UM4A( JJ)^RHS( J) 

20 JJ=JJ+1 

DIVISION OF XI BY TRIANGULAR MATRIX 

21 SUM^O.DO 
II^IEND 

DO 2^ I=ISTA,IEND 

RHSC II )=(RHS(II)-SUM)/A( JJ) 

IF(II-ISTA)25,25,22 

22 KK=JJ-1 
SUM^O . DO 

DO 23 J=II,IEND 
SUM = SUM + A(KK)>«RHS( J) 

23 KK=KK+J 
JJ=JJ-II 

2A II^^II-l 

25 IF(IDEF)26,30,26 

26 G0T0(27 , 11 »8) , ISW 

PERFORM STEP (5) 

27 ISW^2 
GOTO 8 

PERFORM STEP (6) 

28 ISTA:::! 
lEND-IRANK 
JJ = 1 
lSl*Jn2 
GOTO 17 

PERFORM STEP (8) 

29 ISUJ = 3 
GOTO 16 

REINTERCHANGE CALCULATED SOLUTION 

30 II“N 
J J--1 
GOTO ^ 

ERROR RETURN IN CASE OF ZERO DIVISOR 

31 IER-1 

32 RETURN 

ERROR RETURN IN CASE OF ILLEGAL DIMENSION 

33 IER--1 
RETURN 
END 


SUBROUTINE TOPLT ( TP , NT , ALOGS , NQ , GS ) 

THIS PROGRAM SPLITS THE XF VECTOR INTO TEMPERATURE AND 
STRESS COMPONENTS FOR PLOTTING AND CURVE FITTING 

DOUBLE PRECISION XF , TP , A L OGS , GS , X , Y , XX , YY , P 
DOUBLE PRECISION XF L T T , T FCN , X^ L TR i , XP L T G , YP L TG 

DIMENSION TP( 32) , ALOGSC 16 ) ,P( 7 KX( 32 ) , Y( AS ) , CTSOT ( 32 ) ,XX( 16 ) , YY( A8) ,GS( 16 ) 
COMMON /ACOMN/ QXPL 1 T( 32) , OXIURT ( 32 ) , QYPL T T ( 32 ) , QXFL TGC 16 ) , QYPLTGC 16 ) 
COMMON /BCOMN/ XP L TT ( 32 ) , T F CN ( 52 ) , XP L T RT ( 32 ) , XP I T G ( 1 6 ) , YP L T G ( 1 6 ) 

COMMON /CCOMN/ XF(A8),NP 
COMMON /DCOMN/ D,E,F,G.H,AA 
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C TEMPERATURE COMPONENT 
C 

5 DO 10 1=1,16 
X(I)=0. 0 
Y(I)=0,0 
10 GS(I)=0.0 

K = 1 

NTfU □ NT-1 
DO 25 I=1,NTM1 
IF (K.EQ.NP) GO TO 20 
15 Y(K) = XFCI) 

XCK)=TP(K) 

K = K+1 
GO TO 25 
20 Y(K)=0.0 

X(K)=TP(K) 

K = K + 1 

IF (I.EQ.NT) GO TO 25 
GO TO 15 
25 CONTINUE 

DO 30 1=1, NT 
XPLTT(I)=X(I) 

XPLTRTC I )=1000 ./(X(I)+A60 . ) 

TFCN(I)=Y(I) 

QXPLTT(I) ° XPLTT(I) 

QYPLTT(I) □ TFCN(I) 

QXPLRT(I) = XPLTRT(I) 

30 CONTINUE 
NB = 2 

CALL CVFITF ( X , Y , NB , NT , P , 0 ) 

D=P(1) 

E=P(2) 

F=P(3) 

C 

C STRESS COMPONENT 
C 

35 NTM2P=NT 

NQMQ=NTM2P+NQ-1 
K = 0 

DO AO I=NTM2P,NQM2 
K = K + 1 

GS(K)=-XF(I) 

AO CONTINUE 
RETURN 
END 

C 

c 

SUBROUTINE CVFITF (XP , YP , NB IG , NSMAL L , P , IPL ) 

C 

C THIS IS A GENERAL CURVE FITTING SUBROUTINE 

DIMENSION LAB(l), KKK(9), L ABH ( 9 ) , L ABX ( 1 ) , L ABY( 1 ) 

DIMENSION X(IOO) ,Y(100),YCALC(100),DELY(100), ERRATA(100),D(8) 

DIMENSION A(S,8) ,P(7) ,CC(7),XP(100) ,YP(100),TITLE(6) ,FMTW(7) ,PK(1),KPL(1A) 
DATA NINC,XNINC/100,99./ 

DATA LABX /’ X ’/ 

DATA LABY /’ F(X) ’ / 

DATA LAB(3)/AH / 

DATA KKK/3,999,2,0,0,0,0,-1,0/ 

C 

C 

C NOTA BENE CVFITF 
C 

C XP ARRAY CANNOT BE LARGER THAN 100 

C YP ARRAY CANNOT BE LARGER THAN 100 

C I.E. NSMALL CANNOT EXCEED 100 

C NBIG CANNOT BE GREATER THAN 6 
C 
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C P □ OUTPUT ARRAY OF COEFFICIENTS WHERE PCD IS THE CONSTANT 

C TERM, P(2) IS THE COEFFICIENT OF THE FIRST DEGREE TERM AND 

C P(3) IS THE COEFICIENT OF THE SECOND DEGREE TERM, ETC. 

C XP = INPUT ARRAY INDEPENDENT VARIABLE 

C YP = INPUT ARRAY DEPENDENT VARIABLE 

C NBIG = DEGREE OF POLYNOMIAL TO BE ESTIMATED 

C NSMALL = NUMBER OF OBSERVATIONS (NO. PAIRS OF DATA PROVIDED) 

C IPL = CVFIT OUTPUT CONTROL 
C 0 FOR NO OUTPUT 

C 1 FOR DATA OUTPUT 

C 2 FOR DATA AND PLOT OUTPUT 

C 

C DETNUM= VARIATION DUE TO REGRESSION 
C DETDEN= TOTAL VARIATION 
C DEV= VARIANCE 

C DVTN= STD. DEVIATION 
C DETRM= COEFFICIENT OF DETERMINATION 
C CORRL= CORRELATION COEFFICIENT 
C ERRATA= RELATIVE ERROR 

C DELY= RESIDUAL ERROR 

C ERS= RELATIVE VARIANCE 

C ERA= RELATIVE STD. DEVIATION 

C 

DATA PK(1) ,KPL(1) ,KPL(2)/1. 0,6^,2/ 

DOUBLE PRECISION A , D , CC , XMAXX , XP , YP , X , P , Y 

1 IF (NSMALL. GT. 100) GO TO 65 
KN-NBIG 

2 IF (KN.GT.6) GO TO 75 
NBIG ° NBIG+1 

4 SN=NSMALL 
DCNBIG+1)=0. 

3 DO 5 1=1, NSMALL 
X( I ) =XP( I ) 

Y(I)=YP(I) 

5 D(NBIG+1)=D(NBIG+1)+Y(I)^)(2 
XMAXX = X( 1 ) 

DO 15 1=2, NSMALL 
IF (XMAXX-DABS(X( I ) ) ) 10,15,15 
10 XMAXX=DABS(X( I ) ) 

15 CONTINUE 
C 

C NORMALIZE THE INDEPENDENT VARIABLE TO THE INTERVAL -l.,+l. 

C 

DO 20 1=1, NSMALL 

20 X( I) =X( I )/XMAXX 
C 

C SET UP THE MATRIX OF COEFFICIENTS 
C 

DO 25 1=1, NBIG 
DO 25 J=1,NBIG 

21 A(I,J)=0. 

DO 25 K=l, NSMALL 

22 FI=FUNC(I,X(K)) 

23 FJ = FUNC( J ,X(K) ) 

24 A( I , J ) =A( I , J )+FI«FJ 
25 CONTINUE 

C 

DO 30 1=1, NBIG 
D(I)=0. 

DO 30 K=l, NSMALL 
FI=FUNC(I ,X(K) ) 

D( I ) =D( I)+Y(K)^FI 
30 CONTINUE 
C 

CALL DTRIA ( A , NB I G , D , DET ) 

CALL DSOLVE (A,NBIG,CC) 

XXAXX=1 ./XMAXX 
DETNUn=0 . 

DO 35 J=1,NBIG 
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P( J)=CC( J) 

DETNUM=DETKUM+CC( J)XD( J) 

35 P( J)=P( J )xfUNC( J ,XXAXX) 

IF (IPL.EQ.O) RETURN 
C 

C COMPUTE Y“CALC,DELTA-Y, AND STANDARD DEVIATION 
C 

DETNUM=DETNUM-SUMYSQ 

DETDEN=D(NBIG+1)-SUMYSQ 

DETRM=DETNUM/DETDEN 

CORRL=SQRT(ABS(DETRM)) 

CN=NBIG 
ENDIV=5N~CN-1 . 

POLREG=DETNUM/ENDIV 
TOTREG=DETDEN/ENDIV 
ERS=0 . 

DEV=0. 

DO 45 I=1,NSMALL 
X(I)=XP(I) 

KDNSMALL+I 
Y(K)=0 . 

DO 40 J=1,NBIG 
FJ = FUNC(J,X(D) 

40 Y(K)=Y(K)+P( J)^FJ 

C 

DELY(I)=Y(I)-Y(K) 

ERRATAC I )=DELY( I )/Y(K) 

ERS=ERS+ERRATA(I)^^2 
45 DEV=DEV+DELY( I)^DELY(I) 

DEV=DEV/ENDIV 

IF (ENDIV.GT.0.0) GO TO 50 
DVTN=1 .000 
GO TO 55 

50 DVTN=SQRT(DEV) 

55 ERS=^ERS/SN 

ERA=SQRT(ERS) 


C 

C 


c 

c 

c 

c 

c 


60 

61 


C 


65 

75 


WRITE (7,85) LABX,LABY 
DO 60 I=1,NSMALL 
K=N5MALL+I 

WRITE (7,90) X(I),Y(I),Y(K),DELY(I),ERRATA(I) 

WRITE (7,95) KN 

WRITE (7,100) (P(I),I=1,NBIG) 

WRITE (7,105) DEV,DVTN,DETRM,CORRL,ERS,ERA 
WRITE (7,110) DET 

IF (IPL.LT.2) RETURN 


RETURN 

WRITE (6,115) N5MALL 
WRITE (6,125) NBIG 
STOP 


C 

85 

90 

95 

100 

105 


110 

115 

120 

125 


FORMAT (1X,3A4,12X,2A4,10X,12HCALC FUNC. ,6X,10H DEVI ATION , 8X, 1~ 
15H RELATIVE ERROR) 

FORMAT (5G18.8) 

FORMAT (41H THE REGRESSION EQUATION FOR THE ABOVE IS/43H Y = A(0)- 
1 + SUM OF (( A(J) ^ X^^^J )), J = 1,I1) 

FORMAT (7G18.8) 

FORMAT (1X,2X,14H THE VARI ANCE=G1 5 . 7 , 2 OH STANDARD DEVI ATI0N=G15 . 7- 
1/3X, 14HDETERMINATI0N=G15.7,8X,12HC0RRELATI0N=G15.7/1X,2X,14H PCT - 
2VARIANCE = G15.7,20H STD. PCT DEVI ATI0N = G15 . 7 ) 

FORMAT (13H0DETERMINANT=G14.6//) 

FORMAT (43HOONLY 100 DATA POINTS ALLOWED, YOUR NSMAL L =I 13 , IH . ) 
FORMAT (6H0NBIG=2A4,2X,7HNSMALL=2A4) 

FORMAT (1X,2X,20H^^POLYNOMIAL DEGREE=I5 , 1 IH IS TOO BIG) 

END 
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c 

c 


SUBROUTINE DTRIA (B,NR,C,DET) 

C 

c 

DIMENSION A(8,8), C(8), B(8,8) 
DOUBLE PRECISION A,B,C 
KM1=NR 
K=NR+1 

DO 5 1=1, KMl 
DO 5 J=1,KM1 
5 A(I, J)=B(I, J) 

DO 10 1=1, KMl 
10 A(1,K)=C(I) 

DO 15 N=2,K 

15 A(1,N)=A(1,N)/A(1,1) 

DO AO J=2,K 

M = 0 

L=J-1 

DO AO 1=2, KMl 
M=M+1 

IF (M-L) 25,25,20 
20 M = L 

25 DO 30 N=1,M 

30 A(I,J)=A(I,J)-A(N,J)XA(I,N) 

IF (I-J) 35, AO, AO 
35 Ad, J)=A(I,J)/A(I,I) 

AO CONTINUE 

DET=1 . 0 
DO A5 1=1, KMl 
DET=DETXA(I,I) 

DO A5 J=1,K 
A5 BCI, J)=A(I, J) 

B(KM1,KM1)=1.0 

RETURN 

END 

C 

C 

SUBROUTINE DSOLVE (A,M,X) 
DIMENSION A(8,8), X(7) 

DOUBLE PRECISION A,X 

MP1=M+1 

MM1=M-1 

DO 10 K=1,MM1 

MMK=M-K 

MMKP1=MMK+1 

1 X(M)=A(M,MP1)/A(M,M) 

SUM=0 . 

DO 5 I=MMKP1,M 
5 SUM=SUM+A(MMK,I)XX(I) 

10 X(MMK)=A(MMK,MP1)-SUM 

RETURN 
END 

C 

c 

FUNCTION FUNC (I,X) 

DOUBLE PRECISION X,Z,ANS 
Z=X 
K = I 

IF (K-8) 15,10,5 
5 WRITE (6,70) K 

STOP 

10 ANS=ALOG10(Z) 

GO TO 65 
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